{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "838ab677",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "96d3aee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def realY(n0,t0,t):\n",
    "    return n0 * np.exp(- t / t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "085993d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "n0 = 1e16 #cm^-3\n",
    "t0 = 1e-6 #s\n",
    "t = np.arange(0,3 * t0,t0 / 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ed6ba8ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x29f111e0700>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEQCAYAAACgBo8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjmElEQVR4nO3deXhV1b3/8fc3EyFkHiAQEhJmGWQwDAIOpQ6oVaxVq7YOSItUadVeb9V7b6ve3t+911qtbZ1rVaytSCta8VZwxgEVgoxhMswJhCRABiABkqzfHzlSjIEc4CT7nJPP63l4yDl77XO+u8t+slh777XNOYeIiIS+CK8LEBGRwFCgi4iECQW6iEiYUKCLiIQJBbqISJhQoIuIhAlPA93MnjGzMjNb5UfbM83sczOrN7PLm23LMbM3zWyNma02s9w2K1pEJEh5PUJ/DpjkZ9utwA3AX1rY9jzwgHPuFGA0UBaI4kREQomnge6c+wDYfeR7ZtbHzOaZ2RIz+9DMBvrabnbOrQAam7UfBEQ5597ytdvrnNvfTocgIhI0vB6ht+Qp4MfOudOAO4DHWmnfH6g0szlmttTMHjCzyDavUkQkyER5XcCRzCweGAf81cy+fLtTK7tFAWcAI2ialnmJpqmZP7ZNlSIiwSmoAp2mfzFUOueGH8c+xcBS59xGADN7FRiLAl1EOpigmnJxzlUDm8zsCgBrMqyV3RYDKWaW4Xs9EVjdhmWKiAQl83K1RTN7ETgbSAd2AvcA7wKPA92BaGCWc+4/zWwU8AqQAtQBpc65wb7PORd4EDBgCTDNOXewfY9GRMRbnga6iIgETlBNuYiIyInz7KRoenq6y83N9errRURC0pIlSyqccxktbfMs0HNzcykoKPDq60VEQpKZbTnaNk25iIiECQW6iEiYUKCLiIQJBbqISJhQoIuIhIlWA721h1D4bs//nZkVmdkKMxsZ+DJFRKQ1/ozQn+PYD6G4AOjn+zONptv2RUSknbUa6C09hKKZycDzrsmnQLKZdQ9Ugc1t2bWP++YWcqihsfXGIiIdSCDm0LOAbUe8Lva99zVmNs3MCsysoLy8/IS+rKhsL89+vJmXlxSf0P4iIuEqEIFuLbzX4opfzrmnnHP5zrn8jIwW71xt1cSBXRmenczv3y3iQH3DCX2GiEg4CkSgFwPZR7zuCWwPwOe2yMz4l/P6U1JZy+zF21rfQUSkgwhEoL8GXOe72mUsUOWc2xGAzz2qCX3TGZ2byu/fLaLukEbpIiLg32WLLwKfAAPMrNjMpprZdDOb7mvyD2AjUAT8Abi5zar9Z0389Lz+lNUc4IVPj7pOjYhIh9LqaovOuatb2e6AWwJWkZ/G9k5jfN80nliwgWvG5BAXE2yPRxURaV8hfafoT88dQMXeg8xcqFG6iEhIB/ppvVL4xoAMnvxgAzV1h7wuR0TEUyEd6NA0Sq/cf4hnPtrsdSkiIp4K+UAf2jOJ8wZ14+mPNlK1X6N0Eem4Qj7QAW4/tz81dfX84cONXpciIuKZsAj0U7onctGp3Xn2403s3nfQ63JERDwRFoEOcPs5/ag91MCTCzZ4XYqIiCfCJtD7dk3g0uFZzPxkM2U1dV6XIyLS7sIm0AF+8s1+HGpwPPaeRuki0vGEVaDnpnfh8pE9+ctnW9lRVet1OSIi7SqsAh3gx9/si8PxyLtFXpciItKuwi7Qe6bE8d1R2cwu2Ma23fu9LkdEpN2EXaADzPhGP8yM373zhdeliIi0m7AM9MykWL4/phdzlpawsXyv1+WIiLSLsAx0gB+d3YfYqAj+9421XpciItIuwjbQMxI6cfM3+vLm6p0s3FDhdTkiIm0ubAMdYOqEPLKSO/Nfr6+hobHF51aLiISNsA702OhI7rxgIKt3VPPykmKvyxERaVNhHegAF5/anZE5yTzw5jr2Hqj3uhwRkTYT9oFuZvz8W4MorznA4+/rZiMRCV9hH+gAI3JSmDy8B3/4cBPFe3SzkYiEpw4R6AB3ThpIhMH989Z5XYqISJvoMIHeI7kz087ozdzl21myZY/X5YiIBFyHCXSAm87qQ9eETvzy9dU06jJGEQkzHSrQu3SK4l/PH8CybZXMXbHd63JERAKqQwU6wHdG9mRIViL3v7GW2oMNXpcjIhIwHS7QIyKMn180iO1VdTz94UavyxERCZgOF+gAY3qnMWlwJo8v2MDOaj1/VETCQ4cMdIC7LxxIfYPj1/N1GaOIhIcOG+i90rpww/hc/vZ5MatKqrwuR0TkpHXYQAeYMbEvKXEx/PL11TinyxhFJLR16EBPjI3m9nP789mm3cwvLPW6HBGRk9KhAx3g6lHZDMxM4L65q7Uao4iENL8C3cwmmdk6Mysys7ta2J5kZnPNbLmZFZrZlMCX2jaiIiP478uGUlpdx4Nv6gSpiISuVgPdzCKBR4ELgEHA1WY2qFmzW4DVzrlhwNnAg2YWE+Ba28zInBS+P6YXMxduZkVxpdfliIicEH9G6KOBIufcRufcQWAWMLlZGwckmJkB8cBuIKTmL/510gDS4jvxb6+spL6h0etyRESOmz+BngVsO+J1se+9Iz0CnAJsB1YCtzrnvpaKZjbNzArMrKC8vPwES24bibHR3HvxYFaVVDPzky1elyMictz8CXRr4b3m1/idDywDegDDgUfMLPFrOzn3lHMu3zmXn5GRcZyltr0Lh2YycWBXHnxzHSWVtV6XIyJyXPwJ9GIg+4jXPWkaiR9pCjDHNSkCNgEDA1Ni+zEz7rtkMM7BPX9fpWvTRSSk+BPoi4F+ZpbnO9F5FfBaszZbgW8CmFk3YAAQkitfZafGcfu5/Xh7TRnzC3d6XY6IiN9aDXTnXD0wA5gPrAFmO+cKzWy6mU33NfslMM7MVgLvAHc65yraqui2NmV8Hqd0T+Te1wqpqTvkdTkiIn4xr6YV8vPzXUFBgSff7Y+lW/dw2eMLuf70XO69ZLDX5YiIAGBmS5xz+S1t6/B3ih7NiJwUrh3bi5mfbGb5tkqvyxERaZUC/RjuOH8AXRM6cfccXZsuIsFPgX4MX16bvnpHNc8t3Ox1OSIix6RAb8WkIZl8c2BXHnxzPcV79ntdjojIUSnQW2Fm3De56aToPX8v1LXpIhK0FOh+6JkSx0/P7c87a8t4Y5XWTReR4KRA99OU8bkMyUrk56+uomLvAa/LERH5GgW6n6IiI3joyuHUHKjnrpdXaupFRIKOAv049O+WwM/OH8Dba3by14Jir8sREfkKBfpxunF8HmN7p3Lf3EK27dZVLyISPBToxykiwvj1FcOIMONfZi+noVFTLyISHBToJ6BnShz3XjKYRZt38/SHIbmopIiEIQX6CbpsZBaTBmfy4JvrWbOj2utyREQU6CfKzPjvy4aS2Dma219axoH6Bq9LEpEOToF+ElK7xPCry4eytrSG37z1hdfliEgHp0A/SRMHduPq0dk8+cEGFm3a7XU5ItKBKdAD4D8uGkR2Shz/8tdl7D1Q73U5ItJBKdADoEunKB66chgle2r55dzVXpcjIh2UAj1A8nNTuemsPrxUsI23Vuvh0iLS/hToAXT7Of05pXsid89ZwS4t4CUi7UyBHkAxURE8/N3hVNfWc+fLK7SAl4i0KwV6gA3ITODuCwfy9poynvpAd5GKSPtRoLeBG8blcuHQTH41f50uZRSRdqNAbwNmxv3fOZWc1Dhm/OVzyms0ny4ibU+B3kYSYqN5/Psjqa47xK2zlmpVRhFpcwr0NjQwM5H/unQoCzfs4uG313tdjoiEOQV6G7v8tJ58Nz+b379bxHvryrwuR0TCmAK9Hdw3eTCndE/k9peWUbxHTzkSkbahQG8HsdGRPP69kTQ0OG75y1IO1jd6XZKIhCEFejvJTe/CA1cMY/m2Sv77H2u8LkdEwpACvR1NGpLJDybk8dzCzcxdvt3rckQkzCjQ29mdFwzktF4p3PXyCorK9npdjoiEEb8C3cwmmdk6Mysys7uO0uZsM1tmZoVmtiCwZYaP6MgIHr1mJLHRkdz85yXsP6j100UkMFoNdDOLBB4FLgAGAVeb2aBmbZKBx4BLnHODgSsCX2r4yEyK5bdXjeCLsr38+yurtIiXiASEPyP00UCRc26jc+4gMAuY3KzNNcAc59xWAOecLrhuxYR+6fz0nP68srSEJxZoES8ROXn+BHoWsO2I18W+947UH0gxs/fNbImZXdfSB5nZNDMrMLOC8vLyE6s4jMyY2JdLhvXg/nlrmbdqh9fliEiI8yfQrYX3ms8RRAGnARcB5wM/N7P+X9vJuaecc/nOufyMjIzjLjbcmBm/uvxURuYkc9tLy1hZXOV1SSISwvwJ9GIg+4jXPYHm19wVA/Occ/uccxXAB8CwwJQY3mKjI3ny2nzSunTiB88vprSqzuuSRCRE+RPoi4F+ZpZnZjHAVcBrzdr8HTjDzKLMLA4YA+juGT9lJHTimRtGse9AA1NnLtaVLyJyQloNdOdcPTADmE9TSM92zhWa2XQzm+5rswaYB6wAFgFPO+dWtV3Z4WdAZgK/v2YEa3ZUc9usZTRquV0ROU7m1SVz+fn5rqCgwJPvDmbPfbyJe+eu5qazenP3Bad4XY6IBBkzW+Kcy29pW1R7FyPHdv24XDaU7+PJBRvpkx7PlaOyW99JRATd+h90zIx7Lh7EGf3S+bdXVvLJhl1elyQiIUKBHoSiIiN49HsjyUvvwvQXlrCpYp/XJYlICFCgB6nE2Gj+eP0oIiOMG59bTOX+g16XJCJBToEexHLS4njy2tMo2VPL9BeWcKC+weuSRCSIKdCD3KjcVH51+al8unE3t764jPoGPe1IRFqmQA8Bl47I4hffGsS8wlLunrNS16iLSIt02WKIuHFCHlW1h/jtO1+Q2Dma/7joFMxaWmZHRDoqBXoIue2cflTVHuKPH20iqXM0P/lmP69LEpEgokAPIWbGL741iJq6eh56az2JsVHcMD7P67JEJEgo0ENMRIRx/3eGUlN3iHvnriaxczSXjezpdVkiEgR0UjQERUVG8LurRzCuTxr/+rcVvFlY6nVJIhIEFOghKjY6kqeuy2dIVhIzXlzKwg0VXpckIh5ToIew+E5RzJwyity0OH44s4Bl2yq9LklEPKRAD3HJcTH8aeoY0uI7ccOzi1i/s8brkkTEIwr0MNAtMZYXpo4hJjKCa//4mRbzEumgFOhhIictjj9NHUN9g+PKJz/hC43URTocBXoYGZCZwKxpYwG46qlPWbOj2uOKRKQ9KdDDTL9uCcy+6XRioiK4+g+fsrK4yuuSRKSdKNDDUF56F2bfdDrxnaK45g+fsmTLHq9LEpF2oEAPU9mpccy+6XTS4mO47o+f8dlGPcpOJNwp0MNYj+TOzL7pdLond+b6Zxfx0Re6+UgknCnQw1zXxFhmTRtLbloXbpy5mPfWlnldkoi0EQV6B5Ae34kXfziW/t3imfanAq39IhKmFOgdREqXGP78g7EM7pHEzX/+nNdXbPe6JBEJMAV6B5LUOZoXfjCGkTkp/OTFpby4aKvXJYlIACnQO5j4TlE8d+Mozuyfwd1zVvLr+etwTs8oFQkHCvQOKC4miqevy+eqUdk88l4RP529nIP1jV6XJSInSU8s6qCiIiP4n8uGkp0axwPz11FaVccT155GUudor0sTkROkEXoHZmbc8o2+/Oa7wyjYspsrnljI9spar8sSkROkQBe+PaInM6eMZkdlHd9+7GMKt2v9F5FQpEAXAMb1TedvPxpHhBlXPvEJC9aXe12SiBwnBbocNiAzgVduHk9OWhdufG4xsxdv87okETkOfgW6mU0ys3VmVmRmdx2j3SgzazCzywNXorSnzKRYZt80lnF90vjZyyt46K31uqxRJES0GuhmFgk8ClwADAKuNrNBR2l3PzA/0EVK+0qIjeaZG0ZxxWk9+d07X/CTWcvYf7De67JEpBX+jNBHA0XOuY3OuYPALGByC+1+DLwMaPWnMBAdGcGvLj+Vn00awOsrtnPZYwvZskvPKhUJZv4EehZw5GRqse+9w8wsC/g28MSxPsjMpplZgZkVlJfrpFuwMzNuPrsvz00ZzY6qOi7+/Ue8v06/r0WClT+Bbi2813xS9WHgTudcw7E+yDn3lHMu3zmXn5GR4WeJ4rWz+mcwd8YEeiR3Zspzi3n0vSLNq4sEIX8CvRjIPuJ1T6D5Un35wCwz2wxcDjxmZpcGokAJDjlpccy5eRwXn9qDB+avY/oLS6ipO+R1WSJyBH8CfTHQz8zyzCwGuAp47cgGzrk851yucy4X+Btws3Pu1UAXK96Ki4nit1cN5+ffGsTba8q49NGPKSrb63VZIuLTaqA75+qBGTRdvbIGmO2cKzSz6WY2va0LlOBiZkydkMcLU8dQuf8Qlz76sR6YIRIkzKu50Pz8fFdQUODJd0tgbK+sZfoLS1hRXMWPJ/bltnP6ExnR0ikXEQkUM1vinMtvaZvuFJUT9uVDqK/M78nv3y3ihmcXUVZd53VZIh2WAl1OSmx0JPd/51T+57KhLN68m0m//ZB31uz0uiyRDkmBLifNzLh6dA6v/3gC3RJjmTqzgF/8fRV1h455FauIBJgCXQKmb9cEXr1lHFMn5PH8J1u45JGPWFta7XVZIh2GAl0CqlNUJD//1iBm3jia3fsOcckjHzNz4WbdiCTSDhTo0ibO6p/BvNvOYHyfNO55rZCpMwuo2HvA67JEwpoCXdpMenwnnrlhFPdePIiPiiqY9PCHenCGSBtSoEubMjNuGJ/HazPGk9olmuufWcS9rxVqOV6RNqBAl3YxMDOR12ZM4IZxuTy3cDPnP/wBHxdVeF2WSFhRoEu7iY2O5N5LBvPStLFERUTwvac/486/raCqVot8iQSCAl3a3Zjeabxx6xlMP6sPf/u8mHMfWqD1YEQCQIEunoiNjuSuCwby6s3jSYvvxLQ/LeGWv3xOeY2uhBE5UQp08dTQnkm8NmM8d5zXn7cKd3LubxYw5/NiXbcucgIU6OK56MgIZkzsxz9unUDv9C78dPZypjy3mJLKWq9LEwkpCnQJGn27JvDX6eO45+JBfLZxN+c+tIBH3yviQL3WhBHxhwJdgkpkhDFlfB5v3n4mE/qm88D8dZz3mw+0gqOIHxToEpSyU+N46rp8nr9xNJERxtSZBUx5dhEby/XIO5GjUaBLUDuzfwbzbj2Tf7/wFBZv3sP5D3/A/76xlr0HdKepSHMKdAl6MVER/PDM3rx7x1lMHp7FEws28M0H3+fVpSW6GkbkCAp0CRldE2L59RXDmHPzOLolxnLbS8u44olPWFVS5XVpIkFBgS4hZ2ROCq/ePJ77vzOUTRX7uPiRj7j9pWVs3bXf69JEPGVe/ZM1Pz/fFRQUePLdEj6qag/x+PsbePbjTTQ6xzWjc5gxsR8ZCZ28Lk2kTZjZEudcfovbFOgSDkqr6vjdu1/w0uJtdIqKYOqEPH54Zm8SY6O9Lk0koBTo0mFsqtjHg2+u4/UVO0iOi+aWs/ty7em9iI2O9Lo0kYBQoEuHs6qkil/NX8cH68vpnhTL7ef057KRWURF6rSRhLZjBbr+65awNCQriedvHM1ffjiGromx/OzlFZz/8Ae8urSE+oZGr8sTaRMKdAlr4/qk8+rN43ji+6cRFRHBbS8tY+KDC3hx0VatESNhR1Mu0mE0NjreXrOTR98rYnlxFZmJsUw7szdXj86hc4zm2CU0aA5d5AjOOT78ooJH3iti0abdpHWJYeoZeVw7thcJuipGgpwCXeQoFm/ezSPvFrFgfTmJsVHcMC6XKePzSOkS43VpIi1SoIu0YmVxFY+89wXzC3cSFxPJlfnZXD8ul7z0Ll6XJvIVCnQRP63fWcMTCzYwd/l26hsdEwd0Zcr4PMb3TcPMvC5PRIEucrzKaur486db+fNnW6jYe5D+3eKZMj6Pb4/I0k1K4qmTDnQzmwT8FogEnnbO/W+z7d8D7vS93Av8yDm3/FifqUCXUHCgvoG5y3fwzEebWL2jmpS4aK4encN1p+eSmRTrdXnSAZ1UoJtZJLAeOBcoBhYDVzvnVh/RZhywxjm3x8wuAO51zo051ucq0CWUOOdYtGk3z368mTdXlxJhxgVDu3Pt2F6Myk3RdIy0m2MFepQf+48GipxzG30fNguYDBwOdOfcwiPafwr0PPFyRYKPmTGmdxpjeqexbfd+nv9kM7MWb2Pu8u30yejCVaNyuGxkFmnxWuVRvOPPCP1yYJJz7ge+19cCY5xzM47S/g5g4Jftm22bBkwDyMnJOW3Lli0nWb6Id/YfrOf/Vuxg1uJtLNmyh+hI47zBmVw1KpvxfdKJiNCoXQLvZEfoLf1X2eJvATP7BjAVmNDSdufcU8BT0DTl4sd3iwStuJgorsjP5or8bNbvrGHWom3MWVrM/63YQXZqZ77r29YtUXPt0j78GaGfTtOc+Pm+13cDOOf+p1m7U4FXgAucc+tb+2LNoUs4qjvUwPzCUl5avI2FG3YRYTBxYFeuyM/m7AEZdIrSFTJyck52hL4Y6GdmeUAJcBVwTbMvyAHmANf6E+Yi4So2OpLJw7OYPDyLzRX7eKlgG38tKObtNWUkxkZx0andmTw8i9G5qZqSkYDz97LFC4GHabps8Rnn3P8zs+kAzrknzOxp4DvAl5Pi9Uf7DfIljdClozjU0MhHRRX8fWkJ8wt3UnuogR5JsVwyPItvj8hiQGaC1yVKCNGNRSJBYv/Bet5avZNXlpbw4RcVNDQ6BmYmcOmILC4Z1oMeyZ29LlGCnAJdJAhV7D3A/63YwStLS1i2rRIzGNUrlUlDMpk0JFPhLi1SoIsEuc0V+/j7su38Y+UO1u2sAWBYzyQmDenOBUMyydUiYeKjQBcJIRvL9zKvsJR5q0pZUVwFwMDMhMMj9wHdEnRnagemQBcJUSWVtcxbVcr8VaUs3rIb5yAvvQvnDe7GxAFdOa1Xih583cEo0EXCQFlNHW+t3sm8VaV8smEX9Y2OxNgozuyfwcSBXTmrf4aWHugAFOgiYaam7hAffVHBu2vLeG9dORV7D2AGw3omM3FgVyYO7MrgHomamglDCnSRMNbY6CjcXs27a8t4d10ZK4orcQ66JnTirP4ZTOiXzrg+6WQkaPQeDhToIh1Ixd4DLFhXzrvryvhwfTnVdfVA04nV8X3TGd83jdF5acR38udGcQk2CnSRDqqh0bGqpIqPN1TwcVEFizfv4WB9I1ERxvDsZF/ApzM8O5mYKJ1cDQUKdBEBmhYPW7JlDx8VVbCwqIIVJVU4B3ExkYzMSWFUbiqj8lIYkZ1C5xgtJBaMTnZxLhEJE7HRkYdH5QBV+w/xycZdLNxQwaJNu3n4nfU4B9GRxtCsJEblpTI6N5X8XqkkxUV7XL20RiN0ETmsav8hlmzdzaJNe1i0aRcrS6o41OAwgwHdEhiVm8ppvVIYnp1Mr7Q4XUXjAU25iMgJqTvUwNKtlSzevJvFm3fz+ZY97DvYAEBKXDTDs5MZnp3CiJxkhmUnk9RZo/i2pikXETkhsdGRnN4njdP7pAFNJ1nX76xh2bZKlm7dw7Jtlby/vpwvx4V9MrowPDuF4TnJDOuZRP9uCcRGay6+vWiELiInpabuECuKq74S8hV7DwIQFWH065bAkB6JDMlKYkhWIqd0TyQuRmPJE6UpFxFpN845ivfUsqqkilXbq1hVUs2qkip27WsKeTPokxF/OOQHdU9kQGaCli3wk6ZcRKTdmBnZqXFkp8ZxwdDuQFPI76w+wMqSKlaVVFG4vYpPN+7m1WXbD++XHt+JU7onMKBbAgMyEzileyJ9u8ZryuY4KNBFpM2ZGZlJsWQmxXLuoG6H3y+vOcDa0mrWldawtrSGtaXV/OnTLRyobwQgwiA3vQsDMxPo3y2Bvl3j6ds1nrz0LnrgdgsU6CLimYyETmQkZHBGv4zD7zU0Ojbv2sfaHTWsK61mbWkNq0qqeWNV6eGTrxEGOalx9O36z5Dv2zWePhldSIjtuFfaKNBFJKhERhh9MuLpkxHPRad2P/x+7cEGNlbspahsLxvK9lJU3vTzgvVlHGr457nAbomdyE3rQl56F3LTu5Cb1oXeGV3ISY0L++kbBbqIhITOMZEM7pHE4B5JX3m/vqGRrbv3U+QL+Y3l+9hcsY+3Vu88fCIWmk7G9kjqTG563OHAz06NIzsljuzUzmExslegi0hIi4qMoHdGPL0z4jmv2bbqukNsrtjHpop9bK7Yz+ZdTT+/vmIHVbWHvtI2JS6anNQ4eqbGkeML+pzUprDvntQ5JBYvU6CLSNhKjI3m1J7JnNoz+WvbKvcfZNvuWrbu3s+2Pfub/t69n8KSKt4sLP3KNI5Z0/ryPZI70yO5Mz19f2cd8Xdi5yjPl0JQoItIh5QcF0NyXAxDeyZ9bVtDo6O0uo6tu5rCfntlLSV7atleVUthSRVvFe7kYEPjV/aJ7xTVdCVPYmzLfyfFkhoXQ0RE24W+Al1EpJnICCPLN/I+nbSvbW9sdOzad5CSytrDYV9SWUtpVR2l1XV89EUFZTV1NDa7bzM60uiWGMv1p+fywzN7B7xuBbqIyHGKiDDfJZedGJ6d3GKb+oZGKvYepLS6rinoq2oprT5AaVUtXRPb5q5YBbqISBuIiow4PNVCdvt8Z/CfthUREb8o0EVEwoQCXUQkTCjQRUTChAJdRCRMKNBFRMKEAl1EJEwo0EVEwoRnzxQ1s3Jgywnung5UBLAcL+lYglO4HEu4HAfoWL7UyzmX0dIGzwL9ZJhZwdEekhpqdCzBKVyOJVyOA3Qs/tCUi4hImFCgi4iEiVAN9Ke8LiCAdCzBKVyOJVyOA3QsrQrJOXQREfm6UB2hi4hIMwp0EZEwEdSBbmaTzGydmRWZ2V0tbDcz+51v+wozG+lFnf7w41jONrMqM1vm+/MLL+psjZk9Y2ZlZrbqKNtDqU9aO5ZQ6ZNsM3vPzNaYWaGZ3dpCm5DoFz+PJVT6JdbMFpnZct+x3NdCm8D2i3MuKP8AkcAGoDcQAywHBjVrcyHwBmDAWOAzr+s+iWM5G3jd61r9OJYzgZHAqqNsD4k+8fNYQqVPugMjfT8nAOtD+P8r/hxLqPSLAfG+n6OBz4CxbdkvwTxCHw0UOec2OucOArOAyc3aTAaed00+BZLNrHt7F+oHf44lJDjnPgB2H6NJqPSJP8cSEpxzO5xzn/t+rgHWAFnNmoVEv/h5LCHB97/1Xt/LaN+f5lehBLRfgjnQs4BtR7wu5usd60+bYOBvnaf7/nn2hpkNbp/SAi5U+sRfIdUnZpYLjKBpNHikkOuXYxwLhEi/mFmkmS0DyoC3nHNt2i/B/JBoa+G95r/d/GkTDPyp83Oa1mjYa2YXAq8C/dq6sDYQKn3ij5DqEzOLB14GbnPOVTff3MIuQdsvrRxLyPSLc64BGG5mycArZjbEOXfkOZuA9kswj9CL+eqzsnsC20+gTTBotU7nXPWX/zxzzv0DiDaz9PYrMWBCpU9aFUp9YmbRNAXgn51zc1poEjL90tqxhFK/fMk5Vwm8D0xqtimg/RLMgb4Y6GdmeWYWA1wFvNaszWvAdb4zxWOBKufcjvYu1A+tHouZZZqZ+X4eTVPf7Gr3Sk9eqPRJq0KlT3w1/hFY45x76CjNQqJf/DmWEOqXDN/IHDPrDJwDrG3WLKD9ErRTLs65ejObAcyn6SqRZ5xzhWY23bf9CeAfNJ0lLgL2A1O8qvdY/DyWy4EfmVk9UAtc5XynwYOJmb1I01UG6WZWDNxD08mekOoT8OtYQqJPgPHAtcBK33wtwL8BORBy/eLPsYRKv3QHZppZJE2/dGY7515vywzTrf8iImEimKdcRETkOCjQRUTChAJdRCRMKNBFRMKEAl1EJACslcXeTuDzcszsTd9CZat9d84ekwJdRCQwnuPrNw6djOeBB5xzp9C0HlRZazso0EVEAqClxd7MrI+ZzTOzJWb2oZkN9OezzGwQEOWce8v32Xudc/tb20+BLiLSdp4CfuycOw24A3jMz/36A5VmNsfMlprZA74blI4paO8UFREJZb4FxsYBf/WtVADQybftMuA/W9itxDl3Pk3ZfAZNq01uBV4CbqBpWYSjUqCLiLSNCKDSOTe8+QbfomMtLaL2pWJgqXNuI4CZvUrTAzCOGeiachERaQO+ZX83mdkVcPhxc8P83H0xkGJmGb7XE4HVre2kQBcRCQDfYm+fAAPMrNjMpgLfA6aa2XKgED+fVOZbR/0O4B0zW0nTuul/aLUGLc4lIhIeNEIXEQkTCnQRkTChQBcRCRMKdBGRMKFAFxEJEwp0EZEwoUAXEQkT/x/3pXrUsK5bGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e78f2118",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calNewW(wi,s,h):\n",
    "    return wi + h *s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0ab08cc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(n,t0):\n",
    "    return - n / t0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9f4f59be",
   "metadata": {},
   "outputs": [],
   "source": [
    "def s_Euler(wi,t0):\n",
    "    return f(wi,t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7cc47694",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_Euler = np.zeros(t.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a621c34f",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_Euler[0] = n0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0fe1d2d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = t0 / 20 # 约在10.5最贴合，步长比这个大，点在线下，步长比这个小，点在线上\n",
    "for i in range(1,t.size):\n",
    "    w_Euler[i] = calNewW(w_Euler[i-1],s_Euler(w_Euler[i-1],t0),h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "d00e95c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEFCAYAAAACFke6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAk2ElEQVR4nO3deXycZbn/8c+Vyd6mSZe0TdMlTZOmSdNmmwkgB2X7CagVRPHY41E5ohUUPR4UBVFBVE4V9bghiKCIC4qIIKIHFVQQ0cxka9KkSfc2aUvTfUnX9P79kcBpY9pMsz2zfN+vV1+vzMzzzFwPN71655r7vh5zziEiIrEtwesARERk9CnZi4jEASV7EZE4oGQvIhIHlOxFROJAotcBDGTKlCkuLy/P6zBERKJKbW3tDudc9kCvRWSyz8vLIxQKeR2GiEhUMbONp3tNZRwRkTigZC8iEgeU7EVE4oCSvYhIHFCyFxGJAyO+GsfM8oHbgEzn3Nv6nksAPg9MAELOuR+O9OcCPFHfyd3PtLFlzyFmZKVx82VFXFWROxofJSISVcKa2ZvZ981su5k193v+cjNrM7M1ZnYLgHNunXPuun5vcSWQCxwDOkYi8P6eqO/k1seb6NxzCAd07jnErY838UR952h8nIhIVAm3jPMQcPnJT5iZD7gHuAIoAZaaWclpzi8CXnLO3QTcMLRQz+zuZ9o4dLSHa/6cRM4OA+DQsR7ufqZtND5ORCSqhJXsnXPPA7v6PV0NrOmbyR8FfkbvDH4gHcDuvp97BjrAzJaZWcjMQl1dXeGEdYotew4xZZ/xusYk7nwojSV/S8LX0/u8iEi8G84XtLnA5pMedwC5ZjbZzO4DKszs1r7XHgcuM7NvAc8P9GbOufudc37nnD87e8Ddvmc0IyuNHZmOW9/XTV1hD299IZnbf5iK/2DaPx37RH0n5y9/jrm3PM35y59TqUdEYt5wvqC1AZ5zzrmdwPX9nuwG+tfxR9TNlxVx6+NN7B/Xw71XHuHvJcd5z+9T+OB3jLXpa8m7Iw9fmu/V2v6hY72/YLxS2wf0Za6IxKzhzOw7gFknPZ4JbBleOEN3VUUu/331InKz0jBgeyAZe3oeM96bw+YvbyZUFmLP83t6a/vHTq0kqbYvIrFuODP7IFBoZnOBTuAdwL+NSFRDdFVF7j/Pzl8LU98xlbZlbTS8roFLKnp49HVwOOXUwwaq7Wspp4jEinCXXj4CvAQUmVmHmV3nnDsO3Ag8A7QCjzrnVo5eqEM38ZKJBFYEmHnTTC5sSOSLD6axeK3vlGNmZJ1a29dSThGJJeGuxlnqnMtxziU552Y65x7se/63zrn5zrl5zrkvjm6ow+Mb56PgqwUcfGg2R1LgpsdSWfZUCuO7IS3Jx82XFZ1yvMo9IhJLIqqfvZktAZYUFBSM2mcsefc8nihK5rlPreG1f/axeGM69ukc3lA+45TjTrdkU+UeEYlGEdUbxzn3lHNuWWZm5qh+zlXnzOLOZy/i3IYAOcUTGHfrNpqvaubIliOvHtO/rHO651XuEZFoEFHJfqyNXzSeipcqmPeVeez+w25qSmrY8sAWnHPcfFkRaUmn1vVV7hGRaBXXyR4gITGBWR+bhX+Fn4yKDNrf307jpY1cljnplKWcuVlp/PfVi/6pPHM25R4REa9EVM3eS+kF6ZQ9W8bWB7ay9ua1BEuDVH1xLn+9+SLMN9D+sV4zstLoHCCxD1QGUm1fRLwS9zP7k1mCMWPZDAIrA0y8ZCJrb1pL3WvqONB84LTnhFvuUW1fRLykZD+A1JmplP66lOKfFnN43WFqK2vZ8LkNnDh64p+O7b9z93TlHtX2RcRLKuOchpkxbek0Jl46kTUfXcOGOzbQ9VgXRd8vYkJgwinHDrhztx8t5RQRL0XUzN7MlpjZ/Xv37vU6lFclZydT8pMSSn9dyrHdx6g7t441H19DT/eAnZpPS0s5RcRLEZXsx2qd/VBMWTKF6pXV5Lw/h46vdhBcHGT3n3cPfmIfLeUUES9FVLKPdImZiRTdV0TZc2UANF7USNsH2ji+9/ig54Zb29dSThEZDarZD8HEi3obq224fQObv7aZnU/vZP6985myZMoZzwuntn82SzlB9X0RCY9m9kPkS/cx7+55VL5USdLEJJrf3EzLv7VwtOvosN433HIPqL4vIuFTsh+mCdUTqKqtIu+OPLoe66KmuIaXf/oyzrkhvV+45R5QfV9EwqcyzghISE4g7/Y8st+WzarrVtH6zla2P7KdwnsLSZ2ZetbvF065B7ScU0TCp5n9CBq3cByVL1Yy72vz2P3sboIlQbZ8dwvuxNBm+YPRck4RCZeS/QgznzHrv2YRaA6QEcig/fp2Gi9ppHtN94h/lpZziki4IirZR+KmqqFKy0+j7I9lFD1QxP76/YQWhdj0lU2cOP7PLReGajSWcz5R38n5y59j7i1Pc/7y5zT7F4kRNtQvEkeT3+93oVDI6zBGzJEtR2j/YDs7n9xJhj+DogeLGL94/Jh9/vnLnxtwOWduVhov3nLxq49fKfec/FtAWpLvtF8Qi0hkMbNa55x/oNciamYfq1JmpFD6q1JKfl7C4Y2Hqa2qZf3t6zlxZORm+Weico+IKNmPETNj6tunUt1azdSlU9l450ZClSH2/n30S1bavSsiWno5xpImJ1H8cDFTl06l/QPt1L+mnpkfncncz8/FN843+BsM0Ujv3tVSTpHoopm9RyZfMZlAc4AZN8yg4386CC4KsvvZ8BurjQbdiEUkdinZeyhxQiLz75lP+V/KsUSj8dJG2t7fxrE9xzyJRzdiEYldKuNEgKzXZuFv9LPhcxvY/JWTGqtdeebGaqNBN2IRiU2a2UcIX5qPecvnUfWPKpKmJtF8VTMr/3UlR7cPr7HaaNDOXZHoo2QfYTKqMqgKVjH3C3PZ8cQOaopr2PbjbUNurDYatJRTJPpEVLKPpR20w5GQlMCc2+bgb/CTXpTOqnetoumNTRzedNjr0ADt3BWJRtpBG+Fcj6PzO52su3UdZkb+l/KZcf0MLMG8Dm1Q2rkrMra0gzaKmc+Y+eGZBJoDTDhvAqs/tJqGCxvobh/5xmojTeUekcihZB8l0vLSWPzMYop+UMTBpoMEFwfZ9KWRbaw20rRzVyRyaOllFDEzcq7NYdJlk1j9odWsu2Ud2x/dTtGDRWSUZ3gd3oC0c1ckMmhmH4VSclIofbyUhY8t5EjnEWr9tay7bR09h3sGPzkCaeeuyOhTso9i2W/Nprqlmmn/Po1Nd22itqKWvX+LvpVM2rkrMvpUxolySZOSKH6omGlLp9G2rI36f6kn98Zc5t41l8Tx0TO8I71zF1TyETmZZvYxYtJlkwg0B8j9UC6d3+4kWBpk1+93eR3WiAp35y6o5CPSn5J9DEnMSKTwW4WUP19OQmoCKy5bwar/WMWx3d40Vhtp4db2QSUfkf6U7GNQ1r9k4W/wM/vW2Wz70TaCJUG6Hu/yOqxhC7e2D1rOKdJf9BR15az4Un3k35VP9jXZtL23jZVvXcmUt06h8NuFpExP8Tq8IQuntg9azinSX0TN7NUbZ+RlVGRQWVPJ3LvmsvM3OwmWBNn2w8hqrDYatJxT5FQRleydc08555ZlZmZ6HUpMSUhKYM6tfY3VStJZde0qVly+gkMbYrekoeWcIqdSGSeOjFswjornK9hy7xbW3bKOYGmQ/P/OJ/dDuVHRWO1s6UYsIv8nomb2Mvoswcj9UC6B5gCZ/5LJmo+sof619RxcddDr0DyhG7FIvFCyj1Opc1JZ/LvFLPjhArpbugmVhdh410ZOHIvcxmqjQZ05JV6ojBPHzIzp757e21jtw6tZf9t6un7R1dtYrTIyG6uNtFfKMIOVZ1TukWinZC8kT0tm4aML6fpVF6s/uJra6lpm3zybOZ+dgy/NN/gbRLmR7MzZ/0Ysr5R7XvkcEa+ojCOvyn5LNoGWANPfM51NyzcRKg+x5697vA4rIqjcI9FOyV5OkTQxiQUPLmDxHxbjjjoaLmig/cZ2ju8/7nVontKNWCTaqYwjA5p06ST8TX7Wf3o9nd/sZOdTO5n/3flMvnyy16F5RjdikWimmb2cVuL4RAq/XkjFixX4xvlouqKJ1ve0cmxnbDRWGw3auSuRSsleBpV5Xib+ej9zPj2H7T/dTk1JDdsf2x7zLReGQjt3JVKpjCNhSUhJYO7n55L9tmxWXbeKlmtamPKWKRTeU0hKTvQ2VhsN2rkrkUgzezkr48vGU/n3SvK/lM+u3+0iWBJk6w+2apZ/lrRzV8aakr2ctYTEBGZ/Yjb+Rj/jFo+j7b1trHj9Cg6t14qTcGkpp4y1iEr2anEcXdLnp1P+p3IK7y1k3z/2ESwN0vGNDlyPZvmDGY2lnE/Ud3L+8ueYe8vTnL/8Oc3+5RQWib9++/1+FwqFvA5DzsLhzYdp/0A7u363iwnnTqDowSLGlYzzOqyod/7y5wZcypmblcaLt1z86uP+O3eh9zeF093JS2KTmdU65/wDvRZRM3uJXqmzUln09CKKf1xM9+puQhUhNnx+AyeOxldjtZGmco+MFCV7GTFmxrR3TqO6pZrsq7PZ8NkN1AZq2Rfa53VoUUs7d2WkaOmljLjkqcmUPFLC1KVTab+hnbpz6ph10yzyPpeHLz32G6uNtJHeuQtazhmPNLOXUTPlzVMIrAyQ894cNn9lM6GyEHv+ssfrsGJSuOUe0HLOeKVkL6MqKSuJou8VUfZsGe6Eo+HCBtpvaOf4vvhurDbSwi33gOr78UplHBkTEy+eSKApwPrPrKfj6x3s/M1O5t83n8lvjN/GaiMtnHIPaPduvNLMXsaML91HwVcLqPxbJb5MH01vaqLl31s4uuOo16HFFe3ejU9K9jLmJpwzAX+dnzm3z6Hr0S6CxUFe/tnLarkwRrScMz4p2YsnEpITmHvHXKpqq0idm0rr0laar2rmSOcRr0OLeVrOGZ9UsxdPjV80nsqXKun4RgfrP72empIa5n1lHjnvy8HMvA4vZulGLPFHM3vxnPmMWTfNItAUIKMqg/Zl7TRe0sihtZpBekk3YoktSvYSMdLmpVH2bBnz75/P/tr9BBcF2fy1zWqs5hHdiCW2qIwjEcXMmPH+GUx+w2Tab2hn7cfWsv3n2yl6sIjxpeO9Di/u6EYssUMze4lIKbkplD5ZSvEjxRxed5jaylo2fE6N1SKRlnJGByV7iVhmxrR3TCPQGiD7mmw23LGB2qpa9tWosVok0VLO6KBkLxEveUoyJT8pofSpUo7tPkbdeXWs+dgaerp7Bj9ZRp1uxBIdVLOXqDHlTVPIasli3SfX0fG1DnY8sYOiB4qYeNFEr0OLeyO5lLP/jVheKfe88jkyNBE1s9dtCWUwiRMSmX/vfMr/XI4lGI0XN9L2gTaO71VjtUinco+3IirZO+eecs4ty8zM9DoUiXBZr8vC3+hn1s2z2PrAVmpKatjx1A6vw5Iz0M5db6mMI1HLl+5j3pfnkf32bNre20bzm5uZunQqBd8oIDk72evwZADaueudiJrZiwzFBP8EqkJV5N2ZR9djXdQU1/DyT9VYLVpp5+7oULKXmJCQnEDeZ/Lw1/tJK0ij9Z2tNC1p4vDmw16HJmdJO3dHh8o4ElPGLRxH5YuVdHyrg/W3rSe4MMi8u+eR8/4cLEGN1aKFdu6OPM3sJeaYz5j10b7GatUZtF/fTsPFDXSv7vY6NBlB2rl7dpTsJWal5adR9ocyih4o4kDDAUKLQ2z6yiZOHFfLhVigpZxnR2UciWlmRs51OUy6YhLtH2xn3c3r6Pp5V29jtcVqrBbNXinDDFaeOdulnLFa8lGyl7iQMiOF0l+V0vWLLlbfuJraqlpmf2o2cz41h4QU/YIbrUZjKWes7t7V/+USN8yMqW+fSnVrNVOXTmXjnRsJVYbY+3ft2I5l4ZZ7ILZLPkr2EneSJidR/HAxi367iJ79PdS/pp41/7WGnoNqrBaLwl3KCbG9e1dlHIlbk6+YTKA5wLpb19Hx9Q52PLmDou8VMfESNVaLNeGUeyC2d+9qZi9xLXFCIvPvmU/5X8qxRKPx0kZWvW8Vx/Yc8zo08UAs795VshcBsl7b11jtk7PY9tA2giVBup7o8josGWOxvHtXZRyRPr40H/OWz2PqNVNZdd0qVr5lJdnXZFP4rUKSp6mxWryI1d27mtmL9JNRlUFVsIq5X5jLjid3UFNSw7YfbVNjNXlVNO7eVbIXGUBCUgJzbpuDv8FPelE6q969iqY3NHF4kxqrSXTu3lWyFzmDccXjqHihgoJvFrDnhT0EFwbpvKcTd0Kz/HgWjTdiUc1eZBDmM2Z+eCaTl0ymfVk7q29czfafbafogSLSi9K9Dk88Em03YtHMXiRMaXlpLH5mMUU/KOJg80GCZUE2Lt/IiWNqrCYDi6SlnEr2ImfBzMi5NodAa4DJb5zM+lvXU3dOHfvr93sdmkSgSFrKqTKOyBCkTE+h9JeldP2yi/YPtVMbqGX2J2cz5zNz8KX6Bn8DiRsjvZRzqDSzFxmG7LdmU91SzfR3TWfTXZsIlYfY+6Iaq8nZCXcp53Ao2YsMU9KkJBb8YAGLn1nMicMnqL+gntUfWc3xA8e9Dk2ixNl05hwqJXuRETLp9ZMINAfIvTGXzm93EiwNsuv3u7wOS6LA2XTmHCqLpF2BZrYEWFJQUPD+1atXex2OyJDtfXEvbe9ro3tVN9Ovnc68r84jaVKS12FJjDOzWuecf6DXImpm75x7yjm3LDMz0+tQRIYl8/xMquqrmH3bbLb9aBs1JTV0/VKN1cQ7EZXsRWKJL9VH/hfyqQpVkTIjhZVvW0nz25o5su2I16FJHFKyFxllGeUZVNZUkr88n52/2UmwJMjWh7aqsZqMKSV7kTGQkJjA7E/OJtAYYNzCcbT9RxsrLlvBoQ3Rf7s7iQ5K9iJjKL0onfK/lFN4TyH7XtpHsDRIx7c61FhNRp2SvcgYswQj94O5BFYGyLogizUfWUP9BfUcbD3odWgSw5TsRTySOjuVRb9dxIKHF9C9qptQeYiNd6mxmowOJXsRD5kZ0981neqWaqZcOYX1t62nrrqO/XVqrCYjS8leJAIkT0tm4aMLWfj4Qo5uO0ptdS3rbl1Hz6GewU8WCYOSvUgEyX5LNoGWANOvnc6m5b2N1fb8dY/XYUkMULIXiTBJE5NY8MACFv9hMe6oo+GCBtpvbOf4fjVWk6FTsheJUJMu7W2sNvOjM9nynS0EFwbZ+b87vQ5LopSSvUgE843zUfA/BVS8WIFvvI+mK5pofU8rx3Ye8zo0iTJK9iJRIPO8TPz1fuZ8Zg7bf7qdmpIatv9iu1ouSNiU7EWiREJKAnPvnNvbWG1WCi1vb2Hl1Ss5slWN1WRwSvYiUWZ82Xgq/15J/pfz2fW/u6gprmHr99VYTc5MyV4kCiUkJjD75tn4V/gZXzaetuvaWPH6FRxar8ZqMjAle5Eoll6YTvmfyim8t5B9/+hrrPaNDlyPZvlyKiV7kShnCUbu9X2N1S7MYs1H11D/L/UcbFFjNfk/SvYiMSJ1ViqLfrOI4h8X0726m1BFiA2f38CJo2qsJkr2IjHFzJj2zmlUt1STfXU2Gz67gdpALftC+7wOTTymZC8Sg5KnJlPySAmlT5ZybMcx6s6pY+0n1tLTrcZq8UrJXiSGTXnzFKpbqsm5LofNd28mVBZiz1/2eB2WeEDJXiTGJWYmUnR/EWXPluFOOBoubKD9hnaO71NjtXiiZC8SJyZePJFAU4CZN81ky/19jdWeVmO1eKFkLxJHfOk+Cr5aQOXfKvFl+mh6UxMt72zhaNdRr0OTUaZkLxKHJpwzAX+dn7w78uj6RRfBkiAv/+xltVyIYUr2InEqITmBvNvzqKqrIjU/ldalrTRf2cyRTjVWi0VK9iJxbnzpeCr/Vsm8r85j9x93U1NSw5bvbdEsP8Yo2YsI5jNm3TSLQFOAjKoM2pe103hJI4fWqrFarFCyF5FXpc1Lo+zZMuZ/bz77a/cTXBRk81c3q7FaDFCyF5FTmBkz3jeD6pZqJl46kbUfX0vdeXUcaD7gdWgyDEr2IjKglNwUSp8spfiRYg5vOExtZS3r71ivxmpRSsleRE7LzJj2jmkEWgJkvz2bjZ/bSKgyxL4aNVaLNkr2IjKo5CnJlPy4hEW/WUTP3h7qzqtjzcfWqLFaFFGyF5GwTX7jZAIrA8xYNoOOr3UQXBRk9592ex2WhEHJXkTOSuKERObfO5/yP5djCUbjxY20LWvj+F41VotkSvYiMiRZr8vCv8LPrE/MYuuDW6kpqWHHr3d4HZachpK9iAyZL83HvC/No/IflSRNTqL5ymZalrZwdLsaq0UaJXsRGbYJ/glUharI+3weXY93UVNSw8s/UWO1SKJkLyIjIiE5gbxP5+Gv95NemE7rv7fStKSJw5sPex2aoGQvIiNsXMk4Kv5aQcHXC9jzpz0EFwbpvK8Td0KzfC+NeLI3s3wze9DMHjvpuQvN7AUzu8/MLhzpzxSRyGI+Y+Z/zuxtrFadweobVtNwUQPdq7u9Di1uhZXszez7ZrbdzJr7PX+5mbWZ2RozuwXAObfOOXddv7dwwAEgFegYicBFJPKl5adR9ocyih4s4kDjAUKLQ2z68iZOHFfLhbEW7sz+IeDyk58wMx9wD3AFUAIsNbOS05z/gnPuCuCTwOeGFqqIRCMzI+e9Ob2N1S6byLpPrqPu3DoONKqx2lgKK9k7554HdvV7uhpY0zeTPwr8DLjyNOe/8s/4biBloGPMbJmZhcws1NXVFVbwIhI9UmakUPqrUkp+XsKRzUeo9dey/jPrOXFEs/yxMJyafS6w+aTHHUCumU02s/uACjO7FcDMrjaz7wI/Ar490Js55+53zvmdc/7s7OxhhCUikcrMmPr2qVS3VDN16VQ2fmEjoYoQe1/a63VoMS9xGOfaAM8559xO4Pp+Tz4OPD6MzxKRGJI0OYnih4uZunQq7R9op/78enI/kkv+F/PxjfN5HV5MGs7MvgOYddLjmcCW4YUjIvFk8hV9jdU+OIPOb3QSLA2y64/9K8YyEoaT7INAoZnNNbNk4B3Ar0cmLBGJF4kZicz/9nzKny/HkowV/28Fq65bxbE9x7wOLaaEu/TyEeAloMjMOszsOufcceBG4BmgFXjUObdy9EIVkViWdUEW/kY/s2+ZzbYfbiNYEqTrCS3WGCkWib0r/H6/C4VCXochIh7ZX7eftuvaONBwgOxrsin8ViHJ05K9DivimVmtc84/0GsR1S7BzJaY2f179+qbeZF4llGZQWVNJXO/OJcdT+6gpriGbQ9vU2O1YYioZO+ce8o5tywzM9PrUETEYwlJCcz51Bz8jX7Si9NZ9Z5VNL2hicOb1FhtKCIq2YuI9DduwTgqXqig4JsF7Hmhr7HaPWqsdraU7EUk4lmCMfPDMwk0B5hw3gRW37iahtc10N2mxmrhUrIXkaiRlpfG4mcWU/SDIg42HyRYFmTj8o1qrBYGJXsRiSpmRs61OQRaA0x+42TW37qeunPq2N+w3+vQIpqSvYhEpZTpKZT+spSFjy3kSGdvY7V1t62j53CP16FFJCV7EYlq2W/Nprqlmunvms6muzYRKg+x90Ut3+4vopK91tmLyFAkTUpiwQ8WsPiZxZw4fIL6C+pZ/ZHVHD9w3OvQIkZEJXutsxeR4Zj0+kkEmgPk3phL57c7CS4MsusZNVaDCEv2IiLDlTg+kcJvFlLxQgUJaQmsuHwFrde2cmxXfDdWU7IXkZiUeX4m/gY/sz81m5d//DI1JTV0/TJ+G6sp2YtIzPKl+sj/Yj5VoSpSZqSw8m0rab66mSNbj3gd2phTsheRmJdR3ttYLX95Pjt/u5NgSZCtP9gaV43VlOxFJC4kJCYw+5OzCTQGGFc6jrb3trHishUc2nDI69DGhJK9iMSV9KJ0yv9STuE9hex7aR/B0iAd3+zA9cT2LD+ikr3W2YvIWLAEI/eDuQRWBsi6IIs1/7mG+gvqOdh60OvQRk1EJXutsxeRsZQ6O5VFv13EgocX0N3WTag8xMYvbuTEsdhrrBZRyV5EZKyZGdPfNZ3qlmqmXDWF9Z9eT62/lv21sdVYTcleRARInpbMwp8vZOGvFnKs6xi159Sy9pa19ByKjcZqSvYiIifJviqbwMoA06+dzuYvbSZUFmLP83u8DmvYlOxFRPpJmpjEggcWUPbHMtxxR8PrGmj/UDvH90VvYzUlexGR05h4yUQCTQFmfnQmW+7dQrA0yM7f7fQ6rCFRshcROQPfOB8F/1NAxd8q8GX4aHpDE63vbuXYzuhqrKZkLyIShsxzM/HX+ZnzmTlsf2Q7NcU1bH90e9S0XFCyFxEJU0JKAnPvnEtVbRWpc1Jp+dcWmt/SzJEtkd9YTcleROQsjV88noqXKsj/cj67n9lNTUkNWx+M7MZqEZXs1S5BRKJFQmICs2+ejb/Jz/jy8bS9r43GSxs5tC4yG6tFVLJXuwQRiTbpBemUP1fO/Pvmsz+4n+CiIJu/vjniGqtFVLIXEYlGlmDM+MAMAi0Bsi7KYu1/raXu/DoOroycxmpK9iIiIyR1ZiqLnlpE8U+KObTmEKGKEBs+v4ETR71vrKZkLyIygsyMaf82jerWarLfms2Gz26g1l/LvuA+T+NSshcRGQXJ2cmUPFJC6ZOlHNt5jLpz61j7ibX0dHvTWE3JXkRkFE158xSqW6rJeV8Om+/ubay2+8+7xzwOJXsRkVGWmJlI0XeLKHuuDOccjRc10nZ9G8f3jl1jNSV7EZExMvGiiQRWBJj5sZls/d5WahbWsPPpsWmspmQvIjKGfOk+Cr5SQOVLlSRNTKLpTU20vLOFo11HR/VzlexFRDwwoXoCVbVV5N2RR9cvugiWBHn5Zy+PWssFJXsREY8kJCeQd3seVXVVpOan0rq0lea3NI/K7tvEEX/HYTCzJcCSgoICr0MRERkz40vHU/m3Sjq+2cHRLUcxn434Z1gkdmnz+/0uFAp5HYaISFQxs1rnnH+g11TGERGJA0r2IiJxQMleRCQOKNmLiMQBJXsRkTigZC8iEgeU7EVE4oCSvYhIHIjITVVm1gVsHMZbTAF2jFA4XoqV6wBdS6SKlWuJleuA4V3LHOdc9kAvRGSyHy4zC51uF1k0iZXrAF1LpIqVa4mV64DRuxaVcURE4oCSvYhIHIjVZH+/1wGMkFi5DtC1RKpYuZZYuQ4YpWuJyZq9iIicKlZn9iIichIlexGROBC1yd7MLjezNjNbY2a3DPC6mdk3+15fYWaVXsQZjjCu5UIz22tmDX1/PutFnIMxs++b2XYzaz7N69E0JoNdS7SMySwz+5OZtZrZSjP7zwGOiYpxCfNaomVcUs2sxswa+67lcwMcM7Lj4pyLuj+AD1gL5APJQCNQ0u+YNwC/Aww4F/iH13EP41ouBH7jdaxhXMtrgUqg+TSvR8WYhHkt0TImOUBl388ZQHsU/10J51qiZVwMGN/3cxLwD+Dc0RyXaJ3ZVwNrnHPrnHNHgZ8BV/Y75krgYdfr70CWmeWMdaBhCOdaooJz7nlg1xkOiZYxCedaooJzbqtzrq7v5/1AK5Db77CoGJcwryUq9P23PtD3MKnvT//VMiM6LtGa7HOBzSc97uCfBz2cYyJBuHGe1/cr3+/MbOHYhDbiomVMwhVVY2JmeUAFvbPIk0XduJzhWiBKxsXMfGbWAGwH/uCcG9VxSRzqiR4b6Nbr/f9VDOeYSBBOnHX09rw4YGZvAJ4ACkc7sFEQLWMSjqgaEzMbD/wS+Khzbl//lwc4JWLHZZBriZpxcc71AOVmlgX8ysxKnXMnf0c0ouMSrTP7DmDWSY9nAluGcEwkGDRO59y+V37lc879FkgysyljF+KIiZYxGVQ0jYmZJdGbHH/inHt8gEOiZlwGu5ZoGpdXOOf2AH8GLu/30oiOS7Qm+yBQaGZzzSwZeAfw637H/Bp4d9832ucCe51zW8c60DAMei1mNt3MrO/nanrHbeeYRzp80TImg4qWMemL8UGg1Tn3tdMcFhXjEs61RNG4ZPfN6DGzNOBSYFW/w0Z0XKKyjOOcO25mNwLP0Lua5fvOuZVmdn3f6/cBv6X32+w1QDfwH17FeyZhXsvbgBvM7DhwCHiH6/u6PpKY2SP0roaYYmYdwO30fvEUVWMCYV1LVIwJcD7wLqCprz4M8ClgNkTduIRzLdEyLjnAD83MR+8/SI86534zmjlM7RJEROJAtJZxRETkLCjZi4jEASV7EZE4oGQvIhIHlOxFREaZDdJYbwjvN9vMft/XFK6lb0fxGSnZi4iMvof4501Tw/EwcLdzrpje/lrbBztByV5EZJQN1FjPzOaZ2f+aWa2ZvWBmC8J5LzMrARKdc3/oe+8Dzrnuwc5TshcR8cb9wIedc1XAx4HvhHnefGCPmT1uZvVmdnff5qwzisodtCIi0ayvmdtrgF/0dXcASOl77WrgzgFO63TOXUZv3r6A3q6fm4CfA9fS20ritJTsRUTGXgKwxzlX3v+FvgZvAzWse0UHUO+cWwdgZk/Qe3OTMyZ7lXFERMZYX2vm9WZ2Dbx6C8KyME8PAhPNLLvv8cVAy2AnKdmLiIyyvsZ6LwFFZtZhZtcB7wSuM7NGYCVh3qGurw/+x4FnzayJ3r733xs0BjVCExGJfZrZi4jEASV7EZE4oGQvIhIHlOxFROKAkr2ISBxQshcRiQNK9iIiceD/A7MuCg5i6O1uAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t),'m-')\n",
    "plt.scatter(t,w_Euler)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c58de2ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def s_Trape(wi,t0,h):\n",
    "    s1 = f(wi,t0)\n",
    "    wNext = wi + h * s1\n",
    "    s2 = f(wNext,t0)\n",
    "    return(s1 + s2) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f27effb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_Trape = np.zeros(t.size)\n",
    "w_Trape[0] = n0\n",
    "h = t0 / 10 # 约在10最贴合，步长比这个大，点在线下方，步长比这个小，点在线上方\n",
    "for i in range(1,t.size):\n",
    "    w_Trape[i] = calNewW(w_Trape[i-1],s_Trape(w_Trape[i-1],t0,h),h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "77960556",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEFCAYAAAACFke6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjbklEQVR4nO3deXRU9d3H8fdvJjOZCUvCLgE0yBIIKltArNaqVQErblVBUQuCik99WlvrU9EKat1pbWurtSgoLrVqiyiiohWQVhEChB0TEtYkYNgSlkwyyczv+SMJjWmAQJbZPq9zOCeZLb/rPXy93PvJ5xprLSIiEt0coV6AiIg0PQ17EZEYoGEvIhIDNOxFRGKAhr2ISAyIC/UC6tK+fXubkpIS6mWIiESUFStW7LHWdqjrubAc9ikpKSxfvjzUyxARiSjGmG1He06ncUREYoCGvYhIDNCwFxGJARr2IiIxQMNeRCQGNHoaxxhzOvAAkGitvbbqMQfwa6A1sNxaO6uxfy7AnMx8ps3PoqDIR3KSl3uHp3LVwC5N8aNERCJKvY7sjTEzjTGFxph1tR4fYYzJMsbkGGPuA7DWbrbWTqj1EVcCXYByIK8xFl7bnMx8Js9eS36RDwvkF/mYPHstczLzm+LHiYhElPqexnkFGFHzAWOME3gOGAmkATcYY9KO8v5UYIm19ufAnSe31GObNj8Lnz/AdYtcdN5jAPCVB5g2P6spfpyISESp17C31i4G9tV6eCiQU3Uk7wf+RuURfF3ygP1VXwfqeoEx5nZjzHJjzPLdu3fXZ1nfUlDko/0Bw/dWu3jkFS+jvnThDFQ+LiIS6xpygbYLsKPG93lAF2NMO2PMC8BAY8zkqudmA8ONMX8EFtf1Ydba6dbadGtteocOdf627zElJ3nZk2iZPLGElb0C/PBfbqbO8pB+2HvCnyUiEm0acoHW1PGYtdbuBSbVerAEqH0ev1HdOzyVybPXcrBFgD9fWcZXaRX86JN4/ud5Q25CLikPpeD0OptyCSIiYashR/Z5QLca33cFChq2nJN31cAuPHHNmXRJ8mKAwiFuzLweJN/amR1P72B5/+UULS4K1fJEREKqIUf2GUAvY0x3IB8YA9zYKKs6SVcN7PLfUcvzoeOYjmTdnsWq760i+c5kNtzgYdoXOYpoikjMqG/08k1gCZBqjMkzxkyw1lYAdwHzgY3A29ba9U231JPX5vttGLJmCF1/3pX8vxQQvDyXdiv8imiKSMww1tpQr+G/pKen26aqOB79Pwu44m3ostfBl2kV/PX7ZRxKgC5JXr6476Im+ZkiIs3BGLPCWpte13Nh1WdvjBkFjOrZs2eT/YxlrX2sGAeXL3Fx+VcuztiawOsXl5HRRxFNEYleYdWNY62da629PTExscl+RnKSl4o4mPPdch4a52NP6yD/876H/5ubQFlBWZP9XBGRUAqrYd8c7h2eitdVGcHM62B59OZS/n5xOX22OFiWtoyClwoIx1NbIiINEVancZpDdeqmujCtc1svFz+dytmt2pB9WzbZt2VT+GYhqS+mMr94n4rVRCQqxNwF2mOxQcvOl3aSe28uFf4gf/+unw8H+LFV//7xupw8cc2ZGvgiEpaOdYE25k7jHItxGJJvT2bI+iFkpQS57lMXv3rdQ5fdKlYTkcimYV8HT1cPT11Rwp9HldKhyMHDr3i58t8qVhORyBVz5+zrK7mNl6VpPjaklHDjZ/Fc/YWb9Kw45l4f6pWJiJy4sDqyN8aMMsZMLy4uDvVSjqR2DibAX0aV8fsfltKyzHDnCw5yfpFDoKTOpmYRkbAUVsO+OXL29VW7WG13uhvzQQ+Sb+tM3m/zyDgrg/2L9h/3c0REwoHSOCdh/8L9ZN2WRWluKZ1v78zGsV6mfaliNREJLaVxGlmbCyuL1br9ohsFL+3EXp5Lh+UqVhOR8KVhf5KcCU56TOvBC5OCHIq33P0PD3e8H0+rEkU0RST8KI3TQMta+Vg+rrJYbdSSymK1Ny4uY2lfRTRFJHzoyL6BkpO8BJzw3nnlTBnnozApyKS5Hn75fgKleaWhXp6ICKBh32A1i9UKOlgevamUdy4pJ3Wrg4y0DAr+UoANht9FcBGJLRr2DVQ7opnc1sslT6UxbP1QWg1pRfakbFZ/fzUlOSWhXqqIxLCwil7WuHnJbZs2bQr1chrMWsuumbvIuScHW2Y5PKktj7cvJO+gIpoi0viOFb0Mq2FfLdxz9ieqrKCMxWPX4Fp0mM2nBJg5soy8jlYtmiLSqJSzD7H45HgeuvQgz11RSvsDDh6a5eWqf7ko9ymiKSLNQ9HLZlJQ7CO/L1XFam6u+tLNkKw4Zo7UrRBFpOnpyL6ZJCd5ATjshRcv9/Pba0vx+OGBNzzk/DyHwGEVq4lI09GwbyY1I5oAa3sEeHSSn/Jrk8j7XR4ZZ2aw/zMVq4lI09BpnGZS+9631Wmc4QO7ULS4iKyJWay+eDWdJ3bm65u9PK1iNRFpRErjhImAL8DWh7ey/Tc7KPZaZl1aRmavylM7Su2ISH0ojRMBnF4nPZ7swfN3BjmQEOSnsz3c+V48rQ6rWE1EGk6nccLM8hY+Vt4Cly11ccWXLvpVFat9laZiNRE5eWF1ZB9OtyUMlepitbnfqSxW29k2yB0feLhvTgKl21WsJiInJ6yGfTjdljBUaqZ2dra3PD62lLeGl9N7h4OMfhnkP5+vYjUROWFhNeyl7mK14U+kMWzDUFqf05pNP97EqgtWUZKtYjURqT+lcSKItZZds3aR+7NcAr4A3R/uzsqLHEz7LFsxTRE5ZhpHF2gjiDGGzuM603Z4Wzb9eBOb79vMgVOCOEeUYTv95/63gAa+iHyLTuNEoPjO8Zwx+wz+OjZI4kGYOsvDDxe7cFUopikiddORfQT7tKuPLybADQvcjFriZnBVsVouimmKyLfpyD6CJSd5OeyFl37g5zfXleKugPvf8HDb4gQqDlWEenkiEkY07CNYzZjmutMD/OpWH5+nBzjnK0PGGRns+2RfiFcoIuFCp3EiWO1ytXadvAycnsrAwy3ImpjFmuFrOGXcKfR4pgfzthb+VwmbLuKKxA5FL6NUoDTAtke2sf3p7QTbOHnxAh9LepQfeV7laiLRR0VoMcjpcXL646czOGMwu1wV3PF3Nz9+N57EQwZQakck1oTVsFc3TuNrNbAVU8aW8M75fgbkOnlshpdz18aBhYIipXZEYkVYDXt14zSNTu28zDunnCnjfRS0C3Lbh/Hc83Y8adYb6qWJSDMJq2EvTaM6tbOzneWJsaW8dkkZvQqc/OKPTvL+mKdiNZEYoGEfA2qWq2Hg64tclM3uTtvzk8j5SQ6Z52dy+OvDoV6miDQhpXFimLWWb177hpy7cwgcDpAyNYXMi51MW6BiNZFIpDSO1MkYwym3nMLQjUNpf2V7tjywhYNXZBP3dRmW/xSrzcnMD/VSRaSBNOwFdyc3/d7uxxs3BWl1CKa86uHaz124yhXRFIkW+g1aOeKfXXx8MRHGLHBz+VduBmdXFqvlqFhNJOLpyF6OSE7yUuKBmZf5eXq0j7gAPPCGlzs+T6DioIrVRCKZhr0cUbNYbUNKkF/d6mPB0ArOXlpZrLb3470hXqGInCydxpEjaherte/oZdALqQwqbUnWhCzWjlxLp1s60fOZnszbrmI1kUii6KXUS7AsyLZHt7H9ye0EWjt46UIfX/Yoh8qqHRWriYQBRS+lwRzxDrr/ujuDlw9mp7eC2//h5q45KlYTiRQa9nJCWvZvyZQbS3jrAj9nbXby+EtezlujYjWRcKdz9nLCTmnr5aOzfazsVcH4j+OZ+FE8wzbE8fH1JtRLE5GjCKsje1UcR4bq1M43bS1P3VDKrEvL6LHTwT1/dJL3hzxsIPyuA4nEurAa9qo4jgy1i9WyL3Thn9Oddhe2IefuHDLPy+TwBhWriYQTpXGk0VhrKfxrIZt+uonAwQCn/eo0Vg2PU7GaSDNRGkeahTGGTmM7MXTDUDpc04GtU7ZyaFQ2ro0qVhMJNQ17aXTujm7S3kzjtZuDtCiBKa95GL3QhVvFaiIhozSONJkFyT6WTIDRC92MXOZm4KY4Xh5RRraK1USanY7spclUF6u9PNLPU2N8OCxMftPLpEUJVBxQsZpIc9KwlyZTs1ht42mVxWr/PLuCoRmGjH4Z7J2nYjWR5qJhL02mZkTTAB06eEn/c18GfzkIZ6KTtZevZcNNG/Dv8Yd6qSJRT9FLCYmgP8i2x7ex/fHtxCXGUXRPOx63BRQUK6IpcrIUvZSw43A76P5QdwavGIyvk4MWk3fxw5lBEg8aRTRFmoCGvYRUyzNb8tCNPt68sIx+WyuL1b63Kg6fXxFNkcak6KWEXP4BH3lDYWWvAOM/jmf8/HjO3hjHKyNKQ700kaihI3sJueQkLwC721ieHlPKy8PLSNnl4NGXvex4ZoeK1UQagYa9hFzNiCYGPh9QwSOTyjDntCT3nlxWfmclh9YdCu0iRSKcTuNIyNW+921ykpd7R6dy8YBkCt8qJOd/c1gxaAWnPXAaq0aqWE3kZCh6KWHPv8dPzk9zKPxrIfkdgrw0oowtyUFA974VqUnRS4lo7vZu0t5I49Vbgnh98ODrHsYscKtYTeQE6DSORIyFnX18ORGuX+RmRIaLQZuczBxRRpaK1USOK6yO7HVbQjmW5CQvpfHw6nA/T9zgI2jgvr95uXNhAhXFKlYTOZawGva6LaEcS83UTtapQaaM9/HJORUMWW5YlraMPXP3hHiFIuErrIa9yLHUVaw29Lm+DF46CFc7F+uuWMeGGzfg361iNZHalMaRqBD0B9n+1Ha2/XobztZOej3biyV9ypn2iWKaEjuUxpGo53A7SHkwhfTMdLw9vWwcu5Edo7PwbS/V/W9F0LCXKNOiXwsGfTGIeT8I0nubg8dmeLlgVRzGKqYpsU3RS4k6xmn4+xk+FnU1jPs4nnHz4xm2ofL+twWKaUqM0pG9RKXkJC+7kyzTRpcyc0QZpxY6+PXLXq5f4yVYEQz18kSanYa9RKUjMU0Di/tXcP8EHxtODzLyIweZ52RyaI2K1SS2aNhLVKod02zRzUPKm6mkvZVG6bZSVgxewZapWwiWBZmTmc+5Ty6g+33zOPfJBbqIK1FJ0UuJOeV7y8n5WQ7fvPYNgdPd/Pb8g2zo9J/fwFW5mkQqRS9FanC1c9H31b6c+eGZHNxbxi9ecXPDZ27cVb+LpdSORCMNe4lZ7Ua2Y/L4EhYOrGD4chePzvTSd2vlX4mCIqV2JLpo2EtMa9vJy2uX+nn8Rh9BB/zyLS/jP3Jzerw31EsTaVQa9hLTqlM72d2CPDjex7yz/Zy3No4Hnnexe87uUC9PpNFo2EtMq5naqXDBl1fFUfLaqbTq4mH91etZf/16/N+oWE0in9I4InUIlgfZ8fQOtj6yFWdLJz1/35Ov+lWoWE3CmtI4IifI4XJw2gOnkb4qnYTUBL6+5Wvyr1exmkQuDXuRY2jRtwUD/zWQuaOC9NheWaz2/ZUqVpPIoyI0keMwTsPsNB+fd6ksVrv503iGblSxmkQWHdmL1ENykpc9iZbfXF/KS5eV0XV3ZbHamFVeguUqVpPwp2EvUg81i9X+fWYF90/0sa5XkOHzHaw8eyUHMw+Geokix6RhL1IPtYvVWnb10P2vfej3936UFZSxYsgKNj+wmUBpQMVqEpYUvRRpoPJ95eTek8uuV3YRSHHxzPcOsf4UFatJ81P0UqQJudq66PNyH86afxYHivzcM8vN2H+6iVexmoQRDXuRRtL20rZMHl/CZ4Mr+P6KOB6b4aXfFiegYjUJvbAa9saYUcaY6cXFxaFeishJadfRyxsX+3l8bCn+OLj3bQ8T5rnp4VaxmoRWWA17a+1ca+3tiYmJoV6KyEmpTu3kdA0ydbyP98/x8531cdz/vIvd/1CxmoROWA17kUj3rWK1OFh6RRwlfz2N1qd6WX/tetZdu46yXWWhXqbEIKVxRJpBsCJI3m/z2DJ1C84EJz2e6cHSswIqVpNGpTSOSIg54hyc+stTGbJ6CC36tSBrfBYF12ZRtk3FatI8NOxFmlFCagIDPh/A+1cE6Z7v4NEZXi5eoWI1aXoqQhNpZsZheLevj0VdDOM+dnPTP+M5e2McM1WsJk1IR/YiIZCc5GVfa8sz15Ux/QdldN7r4JFXvNyQqWI1aRoa9iIhULNY7cszKrh/Yglrege59BMHK4eu5OBKFatJ49KwFwmB2sVqrbp46fFGH/rN7od/l58VQ1ewefJmAr5AqJcqUULRS5EwU76/nNx7c9k1Yxfe3l4Kf9mOJwu3K6Ipx6XopUgEcbVx0eelPpz16VkcOlROqwl5XPhOgPgyRTTl5GnYi4Spthe35dd3+JmfXs5FKyuL1c7c7FREU06Khr1IGNtW4uPN7/t57KZSSt1wzzseJs5zU7xTEU05MRr2ImEsOamyLTO3S5Cp43y89x0/wzbE8eTMBArfKSQcr7lJeNKwFwljRyKaQEUcvPvdcp6c4MfTzcOG6zew/pr1lO1UsZocn36DViSMVadups3POpLGuWt0Khc935m83+WxdcpWlvVdRs9nerJ0gIrV5OgUvRSJYCWbSsiamEXx4mI2dg8w49Iy9iRV/p3WvW9jj6KXIlEqoVcCAxYO4L2rgqTkO3hsppdLlsdhgipWk2/TaRyRCGcchjmpPj7vbPjRfDdjP6sqVhupYjX5Dx3Zi0SB6mK1311bxl8uL6XTfgcPv+Jl7AovQb+K1UTDXiQq1CxWW9IvwP0TSliVGuTifzpYMWQFB5YfCPUSJcQ07EWiQO1itdZdvPR6vQ9nvHcG5XvKWXn2SnL/L5dAiYrVYpXSOCJRrqK4gtx7c9n54k68Pb2kvpTKotaHvxXnVEwzOiiNIxLD4hLjSJ2eSv/P+mODllUXrGLNxI3s/can+9/GEA17kRjR5qI2DFk7hH+fZzkv08ljM7z0z6387VzFNKOfhr1IDHEmOJlxbgmP3lSKLx5+9ncPd8yNp1UJFBQpphnNlLMXiTHJSV4242PqOB+XL3ExaomLflsS+PBKi7UWY0yolyhNQEf2IjGmOqYZcMJ755UzdZyPvW0sY/7mYN2V6yjLV7FaNNKwF4kxtWOa9PLQ5r3e9PhtD/b/cz/L0pZR8GIB1lrmZOZz7pML6H7fPM59coEu4kYwRS9F5Ahfro+s27IoWlhEebqXR4YWsaNVxZHnVa4W3hS9FJF68fbw0v+z/vR+sTeBtSX8arqbEcsqi9VAqZ1IpmEvIt9ijCF5YjKTJ/jYkBJgzMJ4HnzdQ5fdlRduldqJTBr2IlKnhG4e/nBNGX8eVUr74spitav+7aJbS2+olyYnQcNeROp07/BUvG4nS9Mqi9WW9Qlw1Rdupr7i4cAyFatFGg17EalTzdTO4QSYd7ODQ39IJqHMsPKcleTck6NitQiiNI6InJCKAxVs/uVmCl4owHO6h9SXUvk8qUTFamFAaRwRaTRxrePo/efeDFg0AOMwrL5oNWsnbGSfitXCmoa9iJyUpO8lkb4mncXnW85d5eTxl7wM2KRitXClYS8iJ83pdfLyOSU8cnMph7xw92wPk96Pp9VhRTTDjYrQRKRBkpO8bMXHQz/ycdlSF1d86aLf1gQ+GqVitXCiI3sRaZCaxWpzv1NZrLa7reX6tx2sHbWW0h2loV6ioGEvIg1Uu1jN9PTQdk5vev6+J0ULi8jol0H+C/nYYPgl/2JJo0cvjTGnAw8Aidbaa6seuwD4NbAe+Ju1dtGxPkPRS5Ho4NvsI+v2LIo+KyLx/ER2/DyJpzZuUUSziTQ4emmMmWmMKTTGrKv1+AhjTJYxJscYcx+AtXaztXZCrY+wwCHAA+Sd+CaISCTynu6l/6f9SZ2RStHKg3iu3cpZ8yswQUU0m1t9T+O8Aoyo+YAxxgk8B4wE0oAbjDFpR3n/v6y1I4FfAg+f3FJFJBIZY+h8a2ee+t8K1nUPMHqRmwdf89Ct0KGIZjOq17C31i4G9tV6eCiQU3Uk7wf+Blx5lPdXFaSyH4iv6zXGmNuNMcuNMct3795dr8WLSOTICvp49uoynruilLYHHEyd5eGaxS4K9yii2RwacoG2C7Cjxvd5QBdjTDtjzAvAQGPMZABjzDXGmL8ArwF/quvDrLXTrbXp1tr0Dh06NGBZIhKOkpO8YCCjb4D7J5awtG8FVyxx89hrCRQvKQ718qJeQ4Z9XeFZa63da62dZK3tYa19ourB2dbaO6y1o493cVZEolN1RBPgsBdevNzPn8b4ae9wkXluJpvu3kTgsIrVmkpDfqkqD+hW4/uuQEHDliMi0ao6dVOzMO2m0amcP70TmydvJv8P+ex9by+9X+zN4nY+Fas1soYM+wyglzGmO5APjAFubJRViUhUumpglzqHdu8/9abj6I5kTchizSVrWN+/gv0XlGE9/0ntVL9fTk59o5dvAkuAVGNMnjFmgrW2ArgLmA9sBN621q5vuqWKSDRL+m4S6avT+fx7lmFrnDw2w8ugbBWrNZZ6Hdlba284yuMfAh826opEJGY5vU5eGVbCghQHt37k5ifveliWWsHrl5RRgFI7DRFWRWjGmFHAqJ49e4Z6KSISIslJXrbh45FbShm5zMWVX7hI26ZitYYKq24ca+1ca+3tiYmJoV6KiIRIzWK1D84pZ8p4H9+0t1z3joO1l62ldLuK1U5GWA17EZHaxWqOHh7azelNz2d7UvSvqmK151SsdqJ0D1oRiRi+rT6yb89m/6f7STyvqlgtS8Vq1XQPWhGJCt4UL2fNP4vUl1PZv/og3uu2MuBjFavVh4a9iEQUYwydx1UWq63uEeC6z91MedXDqd+oWO1YNOxFJCJlB3z86eoy/nRVKW0OGabO8vBDFasdVVhFL0VE6is5yUt+kY/lqQE2nupjzAI3o5a4GZbjovjqYhLPVaqvprA6sjfGjDLGTC8uVgOeiBxb7WK1GT/w8+wNftrFucj8biabfrKJikMVIV5l+FAaR0Qi1pzM/P8qTLu8Vye23L+F/D/lE98tntTpqSzuGBvFasdK42jYi0hUKv6imK8nfI0vy8eSswK8fkEph72Vz3ldTp645syoG/iKXopIzEk8N5H0VeksusAydK2Dx2d4Sc+K3WI1DXsRiVpOj5NZZ5fw8I9K2d/SctccD3e9G0/iIUNBUWyldpTGEZGolpzkZXtVsdqIZS6u/reLvtu8fHx5bBWr6cheRKJadWon6IAPh5Xz4HgfOztafvgPB2uGr8G3NTaO8DXsRSSq1S5Wc/bw0P7d3vR6rhcHlhwg44wM8p7NwwbCL6zSmMIqjVOjz/62TZs2hXo5IhLlSreXkn1HNvs+3kfrc1qTOiOVT0uLIjamqeiliMhRWGv55vVvyLk7h/KDAd4/18/76X4ClcGdiIppKnopInIUxhhOufkUhm4Yyto+Qa5c5GLqLA+n7aocj9ES09SwFxEB3J3c/O6yEp69upRWPsOUVz1ct8iFq5yoiGkqeikiUiU5ycvK3j6+PtXH6IVufrDUzeDsOOZeF+qVNZyO7EVEqlTHNEs88PJIP0+N9hFn4bbpDrJ/nE3FgcgtVtORvYhIleqLsNVpnAMD4on7SS+6vlNG3h/y2Dt3L73/0pt2I9vVWcIWzhdxlcYREamH4q+KyZqQRcmGEvw/aMX9qbvZ4woceT4cUjtK44iINFDisETSV6Zz2oOn4fzoAFNeiGfIRidUHS+He2pHw15EpJ4c8Q66P9Kdh37kY0/rID9+38NP3o0n6WBlv044p3Y07EVETlCwt4dHby7lrQv8nLHFyeMzvJy/Oo7kRG+ol3ZUYTXsdVtCEYkE9w5PJT7eyUdnl/OrW31s7xjk1o/jeejdFvg2h+fRfVgNe2vtXGvt7YmJulGwiISvmuVqu9tYXp9kKLm/I56v/WScmcGO3+8Iu2I1pXFERBpJaV4p2ZOy2TdvH63ObkWfGX341N98xWpK44iINANPVw9nzj2Tvm/0xZfjY9mADBb/bAO79vqwQH6Rj8mz1zInM7/Z16ZhLyLSiIwxdLqxE0M3DmVNmuWKz108NMtL952hLVbTsBcRaQLuDm5+P7KE319TSksfPPiah+sXunCHqFhNdQkiIk0kOcnLql4+sk71cf1CN5ctczNoUxwfXNv8a9GRvYhIE6kuVvPFw6wRfp4c48MBTHzRQdakLCqKm69YTcNeRKSJ1L7/7cH+8bg/6EXXe7qy88WdLOu3jL3z9jbLWhS9FBEJgQPLDpA1IYvD6w7T8caObJrQgqeX5TYooqnopYhImGk9tDWDVwwm5aEUvnm7EMeozXT90o+1TRPR1LAXEQkRh9tBytQUnv1xkMLEIHfOrSxWM8HGj2iGVRrHGDMKGNWzZ89QL0VEpNms8vhYdRNcsiKONocMtuowvDEjmmF1ZK9uHBGJRclJXqwDPhlSwVsXln/r8cYSVsNeRCQWVUc0a/K6nNw7PLXRfkZYncYREYlFte992xSFaRr2IiJh4KqBXZr0/rU6jSMiEgM07EVEYoCGvYhIDNCwFxGJARr2IiIxICyL0Iwxu4FtDfiI9sCeRlpOKEXLdoC2JVxFy7ZEy3ZAw7blNGtth7qeCMth31DGmOVHa36LJNGyHaBtCVfRsi3Rsh3QdNui0zgiIjFAw15EJAZE67CfHuoFNJJo2Q7QtoSraNmWaNkOaKJticpz9iIi8m3RemQvIiI1aNiLiMSAiB32xpgRxpgsY0yOMea+Op43xphnq55fY4wZFIp11kc9tuUCY0yxMWZV1Z8poVjn8RhjZhpjCo0x647yfCTtk+NtS6Tsk27GmIXGmI3GmPXGmJ/W8ZqI2C/13JZI2S8eY8wyY8zqqm15uI7XNO5+sdZG3B/ACeQCpwNuYDWQVus1lwEfAQYYBiwN9bobsC0XAB+Eeq312JbzgUHAuqM8HxH7pJ7bEin7pDMwqOrrVkB2BP9dqc+2RMp+MUDLqq9dwFJgWFPul0g9sh8K5FhrN1tr/cDfgCtrveZK4FVb6SsgyRjTubkXWg/12ZaIYK1dDOw7xksiZZ/UZ1sigrV2p7V2ZdXXB4GNQO3S9IjYL/XclohQ9d/6UNW3rqo/tdMyjbpfInXYdwF21Pg+j//e6fV5TTio7zrPqfon30fGmH7Ns7RGFyn7pL4iap8YY1KAgVQeRdYUcfvlGNsCEbJfjDFOY8wqoBD41FrbpPslUu9UZep4rPb/FevzmnBQn3WupLLz4pAx5jJgDtCrqRfWBCJln9RHRO0TY0xL4B/A3dbaA7WfruMtYbtfjrMtEbNfrLUBYIAxJgl41xhzhrW25jWiRt0vkXpknwd0q/F9V6DgJF4TDo67Tmvtgep/8llrPwRcxpj2zbfERhMp++S4ImmfGGNcVA7HN6y1s+t4ScTsl+NtSyTtl2rW2iJgETCi1lONul8iddhnAL2MMd2NMW5gDPB+rde8D9xSdUV7GFBsrd3Z3Auth+NuizHmFGOMqfp6KJX7bW+zr7ThImWfHFek7JOqNc4ANlprnznKyyJiv9RnWyJov3SoOqLHGOMFLga+rvWyRt0vEXkax1pbYYy5C5hPZZplprV2vTFmUtXzLwAfUnk1OwcoAcaHar3HUs9tuRa40xhTAfiAMbbqcn04Mca8SWUaor0xJg+YSuWFp4jaJ1CvbYmIfQKcC9wMrK06PwxwP3AqRNx+qc+2RMp+6QzMMsY4qfwf0tvW2g+acoapLkFEJAZE6mkcERE5ARr2IiIxQMNeRCQGaNiLiMQADXsRkSZmjlOsdxKfd6ox5pOqUrgNVb9RfEwa9iIiTe8V/vuXphriVWCatbYvlf1ahcd7g4a9iEgTq6tYzxjTwxjzsTFmhTHmX8aYPvX5LGNMGhBnrf206rMPWWtLjvc+DXsRkdCYDvyvtXYw8Avg+Xq+rzdQZIyZbYzJNMZMq/rlrGOKyN+gFRGJZFVlbt8B3qlqdwCIr3ruGuCROt6Wb60dTuXc/i6VrZ/bgbeAcVRWSRyVhr2ISPNzAEXW2gG1n6gqeKursK5aHpBprd0MYIyZQ+XNTY457HUaR0SkmVVVM28xxlwHR25B2L+eb88A2hhjOlR9fxGw4Xhv0rAXEWliVcV6S4BUY0yeMWYCMBaYYIxZDaynnneoq+rB/wXwmTFmLZW99y8edw0qQhMRiX46shcRiQEa9iIiMUDDXkQkBmjYi4jEAA17EZEYoGEvIhIDNOxFRGLA/wNbim5UxQhegAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t),'m-')\n",
    "plt.scatter(t,w_Trape)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5a3aa66a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK4(wi,t0,h):\n",
    "    s1 = f(wi,t0)\n",
    "    s2 = f(wi + s1 * h / 2,t0)\n",
    "    s3 = f(wi + s2 * h / 2,t0)\n",
    "    s4 = f(wi + s3 * h / 2,t0)\n",
    "    return(s1 +  2 * s2 + 2 * s3 + s4) / 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "36b40d35",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_RK4 = np.zeros(t.size)\n",
    "w_RK4[0] = n0\n",
    "h = t0 / 1 # 约在10.1最贴合，步长比这个大，点在线下方，步长比这个小，点在线上方\n",
    "for i in range(1,t.size):\n",
    "    w_RK4[i] = calNewW(w_RK4[i-1],RK4(w_RK4[i-1],t0,h),h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ac97dbbc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEFCAYAAAACFke6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZzUlEQVR4nO3db4wc933f8fdn947iWRJ1ism21pIyqYBRqopo6F4YtUYMprAixhBCglFQUkFbp4QJFaDbPIgQqQkqtEVABSoKWLZamTKvNIuWrMuqDJPQYYMkBfPASXk07UqUSoNm6nCPqXmWSkq0TiTv7tsHu0fuLffu5nZm/8zu5wUQvv3N7MxvPNB3Zmc+8xtFBGZm1tsKne6AmZm1nou9mVkfcLE3M+sDLvZmZn3Axd7MrA8MdLoDjaxcuTLWrl3b6W6YmeXK6dOnfxARqxpN68piv3btWsbGxjrdDTOzXJH0vfmm+TKOmVkfcLE3M+sDLvZmZn3Axd7MrA+42JuZ9YHM0ziSHgJ+HbgvIp6qthWAfwWsAMYi4qtZrxfg6JlxXjpxjktXJnlgeIhnn3iYbRtLrViVmVmuJDqzlzQq6bKkN+vat0g6J+m8pOcAIuJCROyqW8RWoATcBMpZdLze0TPjPP/6G4xfmSSA8SuTPP/6Gxw9M96K1ZmZ5UrSM/sDwJeAg7MNkorAK8DjVAr4KUnHIuKtBt9/GPhGRHxZ0hHgD1P1uoGXTpzjwxvT/MKfDDK5DCbvCj5cFhz7i3Ns3nU3xRVFBu4bYGDFAMUVRQqDvoJlZv0jUbGPiJOS1tY1bwLOR8QFAEmHqZzBNyr2ZeBG9e/pRuuQtBvYDfDggw8m6dYcl65MsuwmfOZPBymG5kz71sFv3TF/YahQOQBUi3/tgWDgvoHE0woDPmiYWfdLc82+BFys+VwGfkrSR4HfBDZKej4i9gKvA1+U9NPAyUYLi4h9wD6AkZGRJb9R5YHhIcavTLLr2Q8YnIKPXBfLb8DHly3ny9s3MnV1iun3ppl6b6ry99Vppt6v/u97lWmT35289ffU1SmYWXy9hY8Ubh8IVgxQvG/xg8ScA8Z9RQbuHUBFLb4yM7MmpSn2japTRMQ7wDN1jR8A9dfxM/XsEw/z/OtvMHlzmpuDcHUwuDFY5OntD3P/xvuXvLyIYOaDmcqB4f1K8b91IHhv7kFi6urUnPkmvz85Z95EB427C40PBCsWOEjUT7un6IOGmTWUptiXgTU1n1cDl9J1p3mzqZus0jiSKN5dpHh3MVW/IoLpH07fOijcOjjUHjhqp9V8vn7p+q0DyPT705Dg907xnuL8vy4W+uVRO+2eIir4oGHWS9IU+1PAeknrgHFgB/B0Jr1q0raNpUTFvZ0RTUkM3DPAwD0D3PXAXU0vJ2YqB41GB4XaXxWNDiDXy9dvXbqavtbwlkldp6F4b3HJl6fqpxXvKSL5oGHWDRIVe0mHgM3ASkll4IWI2C9pD3ACKAKjEXG2ZT3NyGxEc/JmpejNRjSBrs7kqyAG7h1g4N50j0bETDB9bfrOexjz/Mq49feVKa5/7/qteRIfNFLeBB+4b4DCRwo+aJillDSNs3Oe9uPA8Ux71GIvnTh3q9DPmrw5zUsnznV1sc+KCqoU0RUpDxrTUbnB3egeRv0B5OrtS1E337nJ5J9P3po288MkNzRY/Eb3YjfBVwxQGPJBw/pXV45n30qXrkwuqd0aU1EMDg8yODyYajkzUzNMvz995yWpRjfBaw4mNyduMvndyVvTZiYTHDSKpL8JvqJIYbkPGpY/fVfsZyOajdqt/QoDBQr3Fxi8P+VB4+bMvPctFroJfuP/3mDyO5O3ps18uPhBQwNKfxN8RZHi8nQ3/82Wou+KfW1Ec9bQYJFnn3i4g72ytAqDBQo/UmDwR1IeNG7MLJiUmvcm+KXrTL19e1pcXzw6pWVKfRN84L4BCsv8YJ8tru+K/VIimh5Yrf8UlhVYtnIZrEy3nJnrM/Pfw1jgAHL94vU588TNBAeNu9TUQeKOy1MeQqSnKWLJD6u23MjISHT6HbT1qR2o/ALYu32DC761zfSH08lvgs9zf2P66jQxtfh/54XlhdQ3wT2ESGdJOh0RI42m9d2ZfVL9ntqx7lBcXrm2v+yvLGt6GRHBzIczC97DaHTgmLo6VRlCpGZaoqfBhwqpb4IX7/VBI2su9vNwasd6hSSKQ0WKQ0WW/dWUB43qECKNhg6541JVzbTJ70/O+V6Sp8FnhxBJcxPc407d5mI/D6d2zOaqHUIki6fBk6Sm6g8g1y9dv/15KUOILHE02zum9cAQIi7283Bqx6w1ap8Gv6uU8qBxLcH9jAbjUl0v1zwN/n6Cp8GpDCGS9iZ48e7OHTRc7OeR9cBqZpatOU+Dr25+OTF9+6BRe1Cof27jjvsaV6b48Hsf3pov0dPg9UOINLiHsfzjy1n9+RQbNN+qncZJzxFNM5uZmmH62nTjXxmz78+Y5+G+2vmHfnSIn/z2TzbVB6dxWiivA6uZWbYKAwUKw4XUQ4jETGtOwJ1tSmmhiKaZ2VK16pq+i31KjmiaWR642Kc0XxTTEU0z6yYu9ik9+8TDDA3OHb3QEU0z6zaZF3tJD0naL+lIXfvdkk5LejLrdXbSto0l9m7fQGl4CAGl4SGPn2NmXSfpawlHgSeByxHxaE37FuALVF5L+JWIeDEiLgC76os98GvA17LpdnfpxnffmpnVSnpmfwDYUtsgqQi8Avwc8AiwU9Ijjb4s6dPAW8D3m+5pzs1GNMevTBLcjmgePTPe6a6ZWR9IVOwj4iTwbl3zJuB8RFyIiBvAYWDrPIv4GeAx4Gngc5LuWK+k3ZLGJI1NTEwk3oC8cETTzDopzTX7EnCx5nMZKEn6qKRXgY2SngeIiF+PiF8B/hPwWkTc8VxxROyLiJGIGFm1alWKbnUnRzTNrJPSPEHbKPkfEfEO8EyjL0TEgRTryzWPomlmnZTmzL4MrKn5vBq4lK47vcsRTTPrpDRn9qeA9ZLWAePADirX5K0Bv/vWzDopafTyELAZWCmpDLwQEfsl7QFOUIlejkbE2Zb1tAckiWh6YDUza4VExT4ids7Tfhw4nmmP+pzffWtmreDhErqMUztm1gou9l3GA6uZWSu42HcZp3bMrBX8pqou43ffmlkruNh3oaQDq4FjmmaWjIt9jjmmaWZJ+Zp9jnlwNTNLysU+xxzTNLOkXOxzzDFNM0vKxT7HHNM0s6R8gzbHPLiamSXlYp9zHlzNzJLwZZw+4NSOmbnY9wGndszMxb4POLVjZi72fcCpHTNrSbGX9JCk/ZKO1LRtk/SapN+W9LOtWK81tm1jib3bN1AaHkJAaXiIvds3+OasWR9RRCSbURoFngQuR8SjNe1bgC9QeTXhVyLixZppRyLiqbrl3A/864jYNd+6RkZGYmxsbEkbYtlwRNMsvySdjoiRRtOWcmZ/ANhSt+Ai8Arwc8AjwE5JjyyynN+ofse6zGxEc/zKJMHtiObRM+Od7pqZpZS42EfESeDduuZNwPmIuBARN4DDwNZG31fFbwFfj4hvNpi+W9KYpLGJiYnkW2CZcUTTrHelvWZfAi7WfC4DJUkflfQqsFHS89Vpnwc+DTwl6Zn6BUXEvogYiYiRVatWpeyWNcMRTbPelfYJWjVoi4h4B3imrvFl4OWU67MWemB4iPEGhd0RTbP8S3tmXwbW1HxeDVxKuUzrEEc0zXpX2jP7U8B6SeuAcWAH8HTqXllH+P23Zr0rcbGXdAjYDKyUVAZeiIj9kvYAJ6hEL0cj4mxLemptkfT9t45omuVL4mIfETvnaT8OHM+sR9b1PIqmWf54uARbMkc0zfLHxd6WzBFNs/xxsbcl8yiaZvnjYm9L5oimWf74tYS2ZH73rVn+uNhbU/zuW7N88WUcaxmndsy6h4u9tYxTO2bdw8XeWsapHbPu4WJvLePUjln38A1aaxkPrGbWPVzsraU8sJpZd3Cxt45zRNOs9XzN3jrOEU2z1nOxt45zRNOs9VzsreMc0TRrvcyLvaSHJO2XdKSm7W5JX5X0mqRfynqdlm+OaJq1XqJiL2lU0mVJb9a1b5F0TtJ5Sc8BRMSFiNhVt4jtwJGI+Bzw85n03HrGto0l9m7fQGl4CAGl4SH2bt8w78Bqn3zxj1j33O/xyRf/iKNnxtvfYbMcSprGOQB8CTg42yCpCLwCPA6UgVOSjkXEWw2+vxp4o/r3dIPp1uc8sJpZayU6s4+Ik8C7dc2bgPPVM/kbwGFg6zyLKFMp+InXaVbPqR2z5qUpvCXgYs3nMlCS9FFJrwIbJT1fnfY68AuS/h3wO40WJmm3pDFJYxMTEym6Zb3KqR2z5qV5qEoN2iIi3gGeqWv8IfDLCy0sIvYB+wBGRkYiRb+sRz0wPMR4g8Lu1I7Z4tKc2ZeBNTWfVwOX0nXHbH5O7Zg1L82Z/SlgvaR1wDiwA3g6k16ZNeCB1cyal6jYSzoEbAZWSioDL0TEfkl7gBNAERiNiLMt66kZyQdWAw+uZlYrUbGPiJ3ztB8HjmfaI7MMOKZpNpdjkNaTHNM0m8vF3nqSY5pmc7nYW0/y4Gpmc7nYW09yTNNsLr+pynrSUmKaTu1YP3Cxt57lwdXMbvNlHOtrTu1Yv3Cxt77m1I71Cxd762tO7Vi/cLG3vubUjvUL36C1vubB1axfuNhb30s6uJojmpZnLvZmCTiiaXnna/ZmCTiiaXnnYm+WgCOalncu9mYJOKJpedeWYi/pQUnHJI1Keq4d6zTLkiOalndNF/tq4b4s6c269i2Szkk6X1PYfwz4vYj4R8AjKfpr1hHbNpbYu30DpeEhBJSGh9i7fYNvzlpuKCKa+6L0KeAacDAiHq22FYHvAI8DZSovJd8JfB84AgTwHyLi3y+07JGRkRgbG2uqX2ad5oimdYqk0xEx0mha09HLiDgpaW1d8ybgfERcqK74MLAVuEnlJeUnJR0BFiz2ZnnliKZ1q6yv2ZeAizWfy9W23wf+iaRXgf/T6IuSdksakzQ2MTGRcbfM2sMRTetWWT9UpQZtERFvAk8t9MWI2Afsg8plnIz7ZdYWjmhat8r6zL4MrKn5vBq4lPE6zLqWI5rWrbIu9qeA9ZLWSVoG7ACOZbwOs67liKZ1q6Yv40g6BGwGVkoqU7kBu1/SHuAEUARGI+JsJj01ywG/+9a6VdPRy1Zy9NJ6XX1qByq/AJzdtzQWil56uASzDnBqx9rNxd6sA5zasXZzsTfrAKd2rN1c7M06wKkdaze/qcqsA/zuW2s3F3uzDvG7b62dXOzNupgHVrOs+Jq9WRdzRNOy4mJv1sUc0bSsuNibdTFHNC0rLvZmXcwRTcuKb9CadTEPrGZZcbE363JJIppO7dhifBnHrAc4tWOLcbE36wFO7dhiXOzNeoBTO7aYthR7SQVJvynpi5L+YTvWadZPnNqxxTRd7CWNSros6c269i2Szkk6L+m5avNWoATcpPJScjPL0LaNJfZu30BpeAgBpeEhv/XK5mj6tYSSPgVcAw5GxKPVtiLwHeBxKkX9FLAT+Hng/0XElyUdiYinFlq2X0to1lqOafamhV5L2HT0MiJOSlpb17wJOB8RF6orPkzlrP4icKM6zzRm1jGOafanrK/Zl6gU9lnlatvrwBOSvgicbPRFSbsljUkam5iYyLhbZjbLMc3+lPVDVWrQFhHxAbBroS9GxD5gH1Qu42TcLzOrckyzP2V9Zl8G1tR8Xg1cyngdZpaCY5r9KetifwpYL2mdpGXADuBYxuswsxQc0+xPaaKXh4BvAA9LKkvaFRFTwB7gBPA28LWIOJtNV80sC45p9qemo5et5OilWXdwRDNfWhK9NLPe5ohmb/HYOGbWkCOavcXF3swackSzt7jYm1lDjmj2Fhd7M2vIEc3e4hu0ZtaQ33/bW1zszWxefv9t7/BlHDNLxamdfHCxN7NUnNrJBxd7M0vFqZ18cLE3s1Sc2skH36A1s1SWktqxznGxN7PUkqR2wBHNTnKxN7O2cESzs3zN3szawhHNznKxN7O2cESzs1zszawtHNHsrLYVe0l3Szot6cl2rdPMuocjmp3V9A1aSaPAk8DliHi0pn0L8AWgCHwlIl6sTvo14Gsp+mpmOeaB1TorTRrnAPAl4OBsg6Qi8ArwOFAGTkk6BjwAvAUsT7E+M8s5D6zWOU1fxomIk8C7dc2bgPMRcSEibgCHga3AzwCPAU8Dn5N0x3ol7ZY0JmlsYmKi2W6ZWc45tdMaWefsS8DFms9l4KciYg+ApM8CP4iImfovRsQ+YB/AyMhIZNwvM8sJp3ZaI+tirwZttwp3RBzIeH1m1mMeGB5ivEFhd2onnazTOGVgTc3n1cCljNdhZj3MqZ3WyPrM/hSwXtI6YBzYQeU6vZlZIh5YrTXSRC8PAZuBlZLKwAsRsV/SHuAElejlaESczaSnZtY3PLBa9pou9hGxc57248DxpntkZpaAI5pL4+ESzCyXHNFcGhd7M8slRzSXxsXezHLJA6stjYu9meWSI5pL4zdVmVkuOaK5NC72ZpZbjmgm52JvZj3NEc0KX7M3s57miGaFi72Z9TRHNCtc7M2spzmiWeFib2Y9zRHNCt+gNbOettSIZq8md1zszaznLSWi2avJHV/GMTOr6uXkjou9mVlVLyd3XOzNzKp6ObnTlmIvaZuk1yT9tqSfbcc6zcyWqpeTO00Xe0mjki5LerOufYukc5LOS3oOICKORsTngM8Cfy9Vj83MWmTbxhJ7t2+gNDyEgNLwEHu3b8j9zVlIl8Y5AHwJODjbIKkIvAI8DpSBU5KORcRb1Vl+ozrdzKwr9ergak2f2UfESeDduuZNwPmIuBARN4DDwFZV/Bbw9Yj4ZqPlSdotaUzS2MTERLPdMjNrudmI5viVSYLbEc2jZ8Y73bV5ZX3NvgRcrPlcrrZ9Hvg08JSkZxp9MSL2RcRIRIysWrUq426ZmWUnjxHNrB+qUoO2iIiXgZczXpeZWUfkMaKZ9Zl9GVhT83k1cCnjdZiZdVQeI5pZF/tTwHpJ6yQtA3YAxzJeh5lZR+Uxotn0ZRxJh4DNwEpJZeCFiNgvaQ9wAigCoxFxNpOempl1iaUMrtYtqR1FRNtXupiRkZEYGxvrdDfMzFKpH1gNKr8AWpXdl3Q6IkYaTfNwCWZmLdJNqR0XezOzFumm1I6LvZlZi3RTasfF3sysRbopteM3VZmZtchSX4nYSi72ZmYt1C0Dq7nYm5l1WDvefetr9mZmHdaOiKaLvZlZh7Ujoulib2bWYe2IaLrYm5l1WDsimr5Ba2bWYe2IaLrYm5l1gaQRzWb5Mo6ZWR9wsTcz6wMu9mZmfcDF3sysD7jYm5n1ga58LaGkCeB7KRaxEvhBRt3ppF7ZDvC2dKte2ZZe2Q5Ity0fj4hVjSZ0ZbFPS9LYfO9hzJNe2Q7wtnSrXtmWXtkOaN22+DKOmVkfcLE3M+sDvVrs93W6Axnple0Ab0u36pVt6ZXtgBZtS09eszczs7l69czezMxquNibmfWB3BZ7SVsknZN0XtJzDaZL0svV6f9L0ic60c8kEmzLZklXJX2r+u+fd6Kfi5E0KumypDfnmZ6nfbLYtuRln6yR9MeS3pZ0VtI/bTBPLvZLwm3Jy35ZLul/Svp2dVv+RYN5st0vEZG7f0AR+C7wELAM+DbwSN08nwG+Dgh4DPizTvc7xbZsBn63031NsC2fAj4BvDnP9Fzsk4Tbkpd98jHgE9W/7wW+k+P/VpJsS172i4B7qn8PAn8GPNbK/ZLXM/tNwPmIuBARN4DDwNa6ebYCB6PiT4FhSR9rd0cTSLItuRARJ4F3F5glL/skybbkQkT8ZUR8s/r3+8DbQP2g6bnYLwm3JReq/19fq34crP6rT8tkul/yWuxLwMWaz2Xu3OlJ5ukGSfv5t6s/+b4u6W+0p2uZy8s+SSpX+0TSWmAjlbPIWrnbLwtsC+Rkv0gqSvoWcBn4g4ho6X7J65uq1KCt/qiYZJ5ukKSf36Qy5sU1SZ8BjgLrW92xFsjLPkkiV/tE0j3AfwV+JSLeq5/c4Ctdu18W2Zbc7JeImAZ+QtIw8N8kPRoRtfeIMt0veT2zLwNraj6vBi41MU83WLSfEfHe7E++iDgODEpa2b4uZiYv+2RRedonkgapFMf/GBGvN5glN/tlsW3J036ZFRFXgP8BbKmblOl+yWuxPwWsl7RO0jJgB3Csbp5jwD+o3tF+DLgaEX/Z7o4msOi2SPprklT9exOV/fZO23uaXl72yaLysk+qfdwPvB0R/2ae2XKxX5JsS472y6rqGT2ShoBPA/+7brZM90suL+NExJSkPcAJKmmW0Yg4K+mZ6vRXgeNU7mafBz4AfrlT/V1Iwm15CvjHkqaASWBHVG/XdxNJh6ikIVZKKgMvULnxlKt9Aom2JRf7BPgk8PeBN6rXhwH+GfAg5G6/JNmWvOyXjwFflVSkckD6WkT8bitrmIdLMDPrA3m9jGNmZkvgYm9m1gdc7M3M+oCLvZlZH3CxNzNrMS0ysF4Ty3tQ0n+vDgr3VvWJ4gW52JuZtd4B7nxoKo2DwEsR8depjK91ebEvuNibmbVYo4H1JP2opN+XdFrSn0j68STLkvQIMBARf1Bd9rWI+GCx77nYm5l1xj7g8xHxt4BfBf5twu/9GHBF0uuSzkh6qfpw1oJy+QStmVmeVQdz+zvAf6mO7gBwV3XaduBfNvjaeEQ8QaVu/zSVUT//AvjPwGepDCUxLxd7M7P2KwBXIuIn6idUB3hrNGDdrDJwJiIuAEg6SuXlJgsWe1/GMTNrs+rQzH8u6Rfh1isI/2bCr58C7pe0qvr57wJvLfYlF3szsxarDqz3DeBhSWVJu4BfAnZJ+jZwloRvqKuOg/+rwB9KeoPKuPevLdoHD4RmZtb7fGZvZtYHXOzNzPqAi72ZWR9wsTcz6wMu9mZmfcDF3sysD7jYm5n1gf8P3LTZI/KAUKMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t),'m-')\n",
    "plt.scatter(t,w_RK4)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "182c31bf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
